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ABSTBACT 


The Langmuir probe technique for measurement of electron concentration 
in the mesosphere is capable of excellent altitude resolution, of order 1 m. 
Measurements from nine daytime rocket flights carrying an electron density 
fine structure experiment frequently show small-scale ionization structures 
in the altitude region 70 to 90 km. The irregularities are believed to be 
the result of turbulent advection of ions and electrons. This work 
describes the fine structure experiment flown by the University of Illinois 
and presents methods of analyzing the collected data. Theories of homo- 
geneous, isotropic turbulence are reviewed. Power spectra of the measured 
irregularities are calculated and compared to spectra predicted by 
turbulence theories. 
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1. INTRODUCTION 


1 . 1 Statement of the Problem 

Electron density irregularities in the mesosphere are often observed by 
both radar and rocket technique. The advection of electrons and ions by 
atmospheric turbulence is one possible cauce of the irregularities. 
Turbulence in the field of a scalar such as electron density may be analyzed 
by decomposing the spatial fluctuations into sinusoidal components. The 
strength of the components as a function of spatial frequency, or wave 
number, is called the scalar Bpectrum. In coherent-scatter radar studies, 
scattered power is measured as a function of aleifude. Irregularities show 
up as isolated regions of high scattered power. Analysis of coherent- 
scatter data gives the value of the scalar spectrum at a particular wave 
number. In rocket studies measuring electron density the irregularities 
show up as small-scale fluctuations on an otherwise smooth electron density 
profile. Analysis of the rocket data can ideally give the scalar spectrum 
at all wave numbers. The measurable wave number range iB limited by system 
noise and performance of the equipment used in the experiment. 

The present study has three main objectives: (1) to review theories of 
isotropic, homogeneous turbulence predicting three-dimensional and one- 
dimensional spectra of turbulence, (2) to describe the analysis of data from 
the electron density fine structure experiment flown by the University of 
Illinois rocket progran, and (3) to compare power spectra of irregularities 
measured by nine rocket flights with theoretical turbulence spectra. 

1 .2 Previous Work 

Rocket studies of mesospheric electron density irregularities have been 
conducted by a number of groups using similar experiments. A rocket 
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equipped with a Langmuir probe passes through an electron density 
irregularity and observes tanporal current fluctuations superposed on a 
background curreut. The tanporal fluctuations are attributed to spatial 
fluctuations in electron density. If the spatial fluctuations are assumed 
to be isotropic, then the temporal and spatial fluctuations are related 
simply by the velocity of the rocket. 

Prakash et al . (1969) describe the results from a number of flights at 
Thumbs, India. Their experiment essentially measures the amplitude of the 
probe current in a frequency range of 70 to 1000 Hz. They observe a number 
of irregularities between 85 and 92 km. The irregularities occur in bands 
with altitude extents of the order of 100 m. Prakash et al. (1969) 
attribute these irregularities to atmospheric turbulence. Power spectra of 
the observed irreguiaritlpi are not calculated. 

Thrane and Grandal (1981.) present fine structure measurements from two 
nighttime launches at the Andoya Rocket Range in Norway. One of the launch 
conditions was the presence of D region irregularities measured with a 
partial reflection radar operated at Tromso, 120 km away from the rocket 
range. The rocket data show ionization fine structure throughout the 
flight. Spectral analysis of the mesospheric irregularities show a power 
spectral index close to -5/3, consistent with the presence of isotropic, 
homogeneous turbulence. 

Royrvik and Smith (1984) compare rocket and radar measurements of a 
single irregularity at 85 km altitude. The rocket launch took place at 
Chilca, Peru. Coincident coherent-scatter radar measurements were made at 
the Jicsmarca radio observatory. Spectral analysis of the rocket data shows 
a power spectral index close to -5/3. The spectrum is extrapolated to the 
wave number measured by the Jicamarca radar, and agreement is found between 
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the rocket and radar results? 

1 .3 University of Illinois Rocket Program 

The l'oc bit program in the Aeronomy Laboratory at the University of 
Illinois conducts three separate experiments based on the Langmuir probe. 
First, the probe current is used to measure the ambient electron density as 
a function of altitude. Calibration between probe current and electron 
density is provided by a radio propagation experiment (Mechtly et al. t 
1967). Second, electron temperature is measured in the E region by sweeping 
the probe potential and analyzing the I-V characteristic of the probe 
(Zimmerman and 8mith , 1980). Third, the ac component of the probe current 
is used to study electron density fine structure. 

The University of Illinois has launched nine daytime flights carrying a 
fine ctructure experiment described by Klaus and Smith (1978). The 
experiment apparatus consists of a Langmuir probe extended from the rocket 
into the ionosphere. As the rocket passes through electron density 
irregularities, temporal fluctuations in the probe current are observed. The 
ac component of the probe current is amplified and telemetered to the 
ground. Digital analysis of the data estimates the power spectra of the 
irregularities . 

1 .4 Outline of the Investigation 

Chapter 2 reviews the theory of the Langmuir probe in the D region. 
Debye length and various plasma parameters in the D region are discussed. 

The values of the parameters in the D region present certain difficulties in 
interpreting the probe current. Some of these problems are overcome by use 
of a dc mode of operation developed by Smith (1969). 

Chapter 3 reviews various theories of isotropic, homogeneous 
turbulence. The basic parameters of fluids and fluid flows are given. The 


spectral approach to describing turbulence and the relationship of three- 
dimensional to one-dimensional spectra are jiscussed. Application of 
dimensional and physical arguments leads to a number of important theories 
predicting turbulence spectra, which are presented. In particular, a scalar 
spectrum model by Hill and Bovhill (1976) suitable for describing 
atmospheric turbulence is discussed. 

Chapter 4 describes the University of Illinois fine structure experi- 
ment and methods of data analysis. Using analog preprocessing, irregulari- 
ties may be quickly identified and distinguished from interference. The 
appropriate data segments are digitized using a microcomputer based 
digitizer. The development and operation of the digitiser are discussed in 
detail. Following transfer of the data to a CDC Cyber mainframe computer, 
the data are Fourier analyzed to calculate power spectra of the irregulari- 
ties. Software needed for the transfer and analysis is described. 

Results of the data analysis are given in Chapter 5. Five of the nine 
flights examined reveal ionospheric irregularities. The power spectra of 
the irregularities are examined and various parameters of the irregularities 
and spectra are tabulated. Comparisons are made to the Hill and Bowhill 
(1976) scalar spectral model and to some spectral predictions by Rastogi and 
Bowhill (1976). 
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2. LANGMUIR PROBE THEORY 


2.1 Introduction 

Rockets using Langmuir probes are capable of measuring electron density 
and other plasma parameters in the ionosphere. In a simple Langmuir probe 
experiment an electrode at a known potential is inserted into a plasma. As 
the potential of the probe is varied the current drawn by the probe is 
monitored. Analysis of this I-V characteristic gives information about 
various properties of the plasma. Due to collisions of particles within a 
few Debye lengths of the probe, exact relations of the plasma parameters to 
probe current do not exist in the D region. Assuming a direct proportion 
between the current and the electron density, a calibration constant between 
the two may be obtained with a radio propagation experiment. When the probe 
is operated in a dc mode high spatial resolution is attainable, allowing 
detailed examination of small-scale ionization structures. Detailed 
discussions of Langmuir probes in the ionosphere are available in Smith 
(1965, 1969) and Cicerone and Bowhill (1967). 

Section 2.2 discusses the D region environment and its relevant 
parameters. Section 2.3 presents the collisionless regime theory and the 
problems it faces in the D region. Section 2.4 describes the dc method of 
operation of the Langmuir probe in the D region. This method gives 
measurements of electron density and fine structure thereof. 

2 .2 The D Region Plasma Environment 

The D region of the ionosphere may be described as a weakly ionized 

plasma consisting of electrons, negative ions, and positive ions. Solar 

radiation is believed to be the chief source of this ionization. Daytime 

2 —3 4 

electron concentrations typically run from 10 cm at 70 km to 10 
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cm -3 at 90 tan; the neutral concentration varies from 10 15 to 10 1 ^ 
cm 3 over the same range. At these densities it is assumed that the 
plasma does not affect the macroscopic dynamics of the neutral atmosphere. 
The free electrons may thus be used as a tracer in detecting macroscopic 
motions of the neutral atmosphere (Thrane and Grandal, 1981). In turn, the 
motion of the neutral atmosphere is assumed not to alter the reaction rates 
of the production and recombination of ions and electrons. 

An ideal plasma consists of a charge-neutral collection of positively 
and negatively charged particles. If the particles have different masses, 
then a gravitational field induces a charge separation, resulting in an 
electric field. Such a separation is present in the earth's iono&phere. 
Lemaire and Scherer (1970) show that in an isothermal ionized atmosphere in 
hydrostatic and diffusive equilibrium the vertical induced field is 
E * g 

where g is the acceleration due to gravity, z is the altitude, and e is the 
electronic charge. y(z) is an effective particle mass: 


y(z) 


h ‘iViVi 
li h\!Vi 


where the sums are over all i species of particles, is the degree of 
ionization (-1 for electrons), m. is the maBS of the i'th species, n. is 
the concentration of the i'th species at altitude z, is the temperature 
of the i'th species, and is Boltzmann's constant. Since the effective 
mass is of the order of the ion mass, the field as given above is quite 
small and has little effect on the probe operation as described in this 


chapter. 

The Debye length is a characteristic length parameterizing the response 
of a plasma to an imposed potential. If a potential is imposed at a point 
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in the plasma the particles rearrange themselves to cancel the resulting 
field. When equilibrium is reached the plasma exhibits a local field in the 
neighborhood of the imposed potential. This field decreases exponentially 
with distance from the imposed potential, in contrast to the inverse square 
decrease of a point charge's field in free space. In the absence of 
negative ions, Poisson's equation in one dimension gives: 

T1 m - <l/e o )e(H . ' V 

dx 

where e is the electron charge, V is the potential, e Q is the permittivity 
of free space, N. is the positive ion density, and N is the electron 
density. The presence of negative ionB would replace K fi by + N_ 
where H, is the negative ion density. Chen (1965) shows thee if the 
positively charged particles are much heavier than the negatively charged 
particles, then 

v - v <r* /h 

o 

where V Q is the imposed potea-l.ial with respect to the plasma, and h is the 
Debye length defined, in rationalized units, as 

h - <e 0 k b T e /N e e 2 ) 1/2 (2.1) 

where T g is the electron temperature, and is the ambient electron 
concentration. At points located a few Debye lengths from the perturbing 
potential the field is small. It is important to note that charge 
neutrality does not hold within a few Deb^s lengths of the imposed 
potential. It is this net charge which shields the rest of the plasma from 
the imposed potential. The charged region is often referred to as the 
plasma sheath. Care must be taken not to visualize the sheath as a sharply 
bounded layer outside of which the field does not exist. The Debye length 
only characterizes the exponential decrease in the field. 
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To calculate the Debye length it is necessary to find the electron 
temperature and the electron density. A number of techniques are available 
for measuring electron density , which may be considered a vellHwown 
paraneter. Given a distribution of velocities for electrons or ions it may 
be possible to apply a Maxwellian distribution, defining an electron or ion 
tanperature. For low collision rates these two temperatures could be quite 
different from each other and that of the neutral atmosphere. Unfortu- 
nately, direct measurements of electron temperature using Langmuir probes 
(Banks and Rockarts, 1973; Zimmerman and Smith, 1980) are obtainable only 
above 100 km. Below this height, conventional Langmuir probe theory fails 
and electron temperature measuronents must be obtained in other manners. 

Due to the high neutral concentration and high collision rates in the D 
region it is assumed that the ion and electrons are in thermal equilibrium 
with the neutral atmosphere (Hill and Bowhill, 1976). Neutral temperatures 
in the D region range from approximately 180 to 220 °K (U.S. Standard 
Atmosphere, 1976). Thus Debye lengths ranging from 1 to 10 cm may be 
expected in the D region. 

2.3 Langmuir Probes in t?ae D Region 

Langmuir probes provide a simple method of measuring charged particle 
densities in the D region. This simplicity, unfortunately, is 
counterbalanced by a rather complex theory. A general probe theory does not 
exist for all atmospheric regions, and useful results may be obtained only 
under certain restrictive assumptions. This section describes the basic 
theory of probes in a plasma and the applicability of the theory to the D 
region. A more complete description of Langmuir probes is available in Chen 
(1965). 

In the basic operation of a Langmuir probe an electrode is inserted 
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into the plasma and a potential imposed. The potential is swept repeatedly 
over a ^ide range, and the induced current is measured as a function of 
potential, the rocket body acting as an electrode for the return current. 
Such an arrangement is termed a bipolar probe. Figure 2.1 shows a typical 
I-V characteristic for a probe. At V , the plasma or space potential, the 
probe potential equals that of the plasma. No sheath forms about the probe 
and the plasma is undisturbed. The observed current is entirely due to the 
thermal migration of charges in the plasma. Since an electron mass is many 
orders of magnitude less than an ion's, the electron thermal velocity is 
much greater than the ion's (for equal ion and electron temperatures). This 
results in a net electron current into the probe. To induce zero current 
the probe potential must be lowered, accelerating positive ions and 
retarding electrons. At V^, the floating potential, the sum of the 
electron and ion currents adds to zero. Alternatively, if the potential is 
raised above V the positive ions are retarded and the electrons are 

o 

accelerated. For potentials much greater than the plasma potential the 
current is due totally to the thermal electrons (region A in Figure 2.1). 

For potentials below V g the current is totally due to positive ions 
(region C in Figure 2.1). The slope of the curve in region B gives 
information about the electron and ion tanperatures. It is important to 
note that an isolated conductor in a plasma assumes the floating potential 
(Chen, 1965). A satellite or rocket in a plasma assumes this potential. 

All measurements of probe potential are thus made with respect to V^. 
Moreover, these measurements must account for the plasma sheath which 
naturally accompanies a conductor not at the plasma potential. 

Quantitative predictions of the I-V behavior of a spherical probe may 
be made for an ideal plasma of electrons and positive ions with Maxwellian 
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distributions. From Smith (1965) the number of electrons collected per unit 
area by a probe at the plasma potential is 

j e - VV 4 (2 - 2) 

where N is the electron concentration, v is the electron mean 
e e 

velocity, and j is the average electron current density. For a 
Maxwellian electron velocity distribution 

v e “ ( 8k b T e /Tnn e ) 1/2 (2.3) 

where T g is the electron temperature, k fa is the Boltzman constant, and 
m 1 b the mass of an electron. A similar analysis applies to positive 
ions. A comparison of the two current densities for equal temperatures 
gives 

i e /ii - ^ i /m e ) lj2 (2.4) 

where j^ is the ion current and m^ is the ion mass. For N0 + , a common 
ion in the D region, this ratio is 235. It is then quite justifiable to 
approximate the measured current as due to electrons only. This 
approximation becomes even better when the probe is above V and positive 
ions are repelled from the probe. If the electron mean free path, L^, is 
much greater than both h and a p , the probe radius, then Smith (1965) shows 
that for electron retarding potentials (V < V fi ) the electron current is 
j - j e expte(V-V s )/k b T e ] (2.5) 

independent of probe shape. A probe operating under these conditions is 
said to be operating in the collisionless regime. For electron accelerating 
potentials on a small spherical probe (h » a ) 

j “ j e ll + e<V-V B )/k b T e l . (2.6) 

Both relations give j proportional to N for constant T . This property 
is essential in justifying the fine structure experiment to be described in 
Section 2.4. From these relations it is theoretically possible to derive 


electron density, electron temperature, and average ion composition. 

There are a number of restrictions on the use of the above relations. 
Many of these restrictions are met in the D region, given due attention to 
the experimental apparatus. First, it must be assumed that the current 
drawn by the probe does not affect the plasma. Given the high electron 
velocities in the D region, any dearth of electrons created locally by the 
probe current is quickly replenished by diffusing electrons. Equation (2.3) 
numerically gives 

v. » 6.21 x 10 3 T lfZ ra/sec 
e e 

or approximately 10 3 m/sec in the D region. Second, photoemission of 
electrons at the probe surface can interfere with current measurenents. 

Smith (1965) finds that photoelectrons in the F region limit electron 
density measurements to altitudes less than 1000 km. It is assumed that 
this effect is absent in the D region since much of the hard solar radiation 
is attenuated by the atmosphere above. Third, it must be assumed that the 
contact potential (the bias appearing between the probe and rocket body when 
the internal circuitry shows no potential) remains constant over the 
altitude range under study. Such a bias might be caused by the dissimilar 
metals in the probe and rocket body (Smith, 1965). Fourth, the telemetry 
system and any radio propagation experiments must be designed such that they 
do not excite any plasma resonances near the probes. 

Two important conditions do not hold in the D region. First, the above 
relations, equations (2.2) through (2.6), assume that the negative current 
is totally due to electrons. However, below 90 km negative ions are present 
in sufficient numbers to disqualify this assumption. Second, the above 
analysis has not taken collisions into account. If the Debye length is 
larger than or of the order of the electron mean free path (h .> l> e ). 
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Chen collisions can occur in the sheath, disrupting the free movement of 
electrons traversing the sheath to the probe* Such a free-fall condition is 
necessary to derive equation (2.6) (Chen, 1965), If the condition L g » 

Sp is violated then collisions can occur in the region close to the probe. 
This region can only receive electrons from that portion of all 4 tt stera- 
dians not subtended by the nearby probe. The anisotropy prevents the 
formation of the Maxwellian velocity distribution necessary to derive 
equation (2.6). Due to the large variation of L , h, and other parameters 
in the D region, it is not possible to derive general relations between 
measured current and electron density. However, Smith (1969) suggests that 
the relationship is a simple proportionality, allowing the current to be 
calibrated by other experiments, thus extending the useful range of the 
Langmuir probe into the D region. 

The measured current is also affected by the rocket's high speed, 
typically on the order of 1 km/sec. This is comparable to the mean positive 
ion velocity (from equation 2.3 NO* at 200° has v + = 664 m/sec) in the 
D region but is far leBS than the mean electron velocity. The effect on 
positive ion current is discussed by Smith (1965) who gives the increase 
expected in ion current as a function of rocket speed divided by v + . For 
example, at v + ® 664 m/sec and rocket velocity of 1 km/sec the NO* ion 
current is 1.6 times the value expected for a stationary rocket. The 
measured electron current is largely unaffected due to the great difference 
between v g and the rocket's velocity. Since the measured ion current 
increases, the total measured current density decreases. However, the ion 
current at zero rocket speed is so small compared to the electron current 
that the effect is negligible. The rocket's motion may also distort the 
current if the probe passes through the rocket's wake. Although the 
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electrons have sufficient velocity to quickly refill the rarefied wake 
region, preservation of the plasma charge neutrality demands that the plasma 
respond to the rarefaction at the slower ion velocity. As the rocket spins 
and precesses, an extended probe may sweep through the trailing rarefied 
wake, modulating the measured current at the spin frequency. Finally, it is 
important to discuss the effect of the supersonic nature of the rocket's 
velocity. Such high velocities are necessary in assuming a frozen 
atmosphere during measurements. At these velocities, however, a shock front 
develops ahead of the probe. This region exhibits temperatures and 
densities above those of the ambient atmosphere. It is found, however, that 
the ratio of the neutral densities n fl h oc j ( V n ®abient de P end ® on *y on t * ie 
Mach number (Thrane and Grandal, 1981) where is the neutral density 

in the shock front and “g^^g^ i- 8 the neutral density in the ambient 
atmosphere. 

2.4 DC Operation of a Rocket- Borne Langmuir Probe 

In the previous section it was found that exact formulae for electron 
density versus current do not exist for D region measurements. Mechtly et 
al. (1967) describe an experiment to calibrate the current. Since 
collisionless regime theory does not hold in the D region, a potential sweep 
of the probe gives questionable results for temperature and other 
quantities. The probe is therefore operated at a fixed potential, 
sacrificing temperature data but gaining altitude resolution. This dc mode 
is particularly useful in discerning any fine structure in electron density. 
These irregularities may be related to irregularities in the neutral 
atmosphere. Such a relationship found by Thrane and Grandal (1981) is given 
here. 

In the dc mode of operation the probe is set at a particular voltage, 
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and the electron current versus time is measured. The main assumption is 
that she current is proportional to the electron density over ranges in 
which the electron temperature is constant (Smith, 1969). The probe 
potential is set high enough above the floating potential that the I-V 
charct eristic is well into the linear part of the curve (Figure 2.1), above 
the plasma potential. The relation of current to electron concentration is 
found by means of a propagation experiment (Hechtly et al., 1967). In this 
experiment a circularly polarized wave at a few megahertz is launched from 
the ground and received at the rocket. The differential absorption and 
Faraday rotation of the wave are monitored as a function of height. These 
quantities can be used to find the electron concentration as a function of 
altitude. The two experiments work in tandem, the uncalibrated probe having 
high spatial resolution and the propagation experiment giving a coarse 
proportionality between concentration and current. 

Amplification of the ac component of the probe current can reveal 
details of the fine structure in electron concentration. Such experiments 
have been carried out by Frakash et al. (1980), Smith and Klaus (1978), and 
Thrane and Grandal (1981). These irregularities can occur over ranges of 
tens of meters to kilometers. To measure them one asumes that the 
fluctuations are statistically stationary, that is, 


N - <N > + AN 
e e e 

where N is the electron concentration and AN is the fluctuation about 
e e 

the mean. If in addition the proportionality between probe current and 
electron density is steady over a particular altitude range, then appropriate 
signal analysis of the current alone yields AN /<N > (■ AI/<I>) (Smith 
and Klaus, 1978) the percentage fluctuation. Due to the proportionality of 
N to I calibration of the current is unnecessary. Thrane and Grandal 
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(1981) extend the proportionality of n B j, ock to and show that 

percentage fluctuations in N g are independent of velocity. 

Further analysis is necessary to relate AH /<N > to An/<n> where n 
is the neutral atmosphere concentration and An is the fluctuation about the 
mean, <n>. The thermodynamic state equation 
p * nkT 

may be differentiated with respect to height in an ideal exponential 
atmosphere, yielding 


i. „ » 1 + JL -dT 

H p H n Tdz (2.7) 

where is the pressure Beale height and is the neutral density 
scale height in the unperturbed atmosphere. The fluctuations are assumed to 
occur over altitude ranges small compared to the scale height. Thrane and 
Grandal (1981) show that for adiabatic displacements of air parcels 

(2 ’ 8> 

where Az is the vertical displacement off equilibrium, and y is the ratio of 
the specific heat of air at constant pressure to that at constant volume. 
Letting H Q be the scale height of the electron concentration, Thrane and 
Grandal (1981) find a relation similar to equation (2.8): 





A _L_ 
V VH 
e 


Thus 

A. ~ 1 

<n> ■ (YH )/H e - 1 <N e > 

giving a relationship between electron density and neutral density 
fluctuations. This analysis assumes that fluctuations are caused only by 
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vertical displacements of air parcels in which the adiabatically shifted 
parcels' densities differ from those of the ambient atmosphere. Ambient 
properties in the horizontal direction are assumed to change slowly with 
respect to the vertical changes. 

2.!) Summary 

Rocket-borne Langmuir probes are capable of measuring a variety of 
plasma properties; electron and ion temperature and electron and ion 
concentrations. Collisionleas regime Langmuir probe theory applies above 90 
km. Below 90 km, in the D region, the Langmuir probe may be used only to 
measure electron density (when calibrated with an independent radio 
propagation experiment) and fluctuations in electron density. Th 'e 
fluctuations in electron concentration may then be related to fluctuations 
in the neutral atmosphere. Homogeneous turbulence is believed to be 
responsible for creating the irregularities. Mechanisms for production and 
a theoretical description of homogeneous turbulence are discussed in the 
next chapter. 


3. THEORIES OF HOMOGENEOUS TURBULENCE 


3,1 Introduction 

Electron density fluctuations as a function of altitude in the 
mesosphere are often observed by both radar and rocket techniques. These 
fluctuations have been observed to occur in layers with a wide range of 
thicknesses ranging from the order of tens of meters (Rastogi and Bowhill, 
1976) to the order of kilometers (Royrvik and Smith, 1984; Rottger et al., 
1979). Isotropic turbulence is generally believed to cause these layers 
(Hocking and Vincent. 1982; Thrane and Grandal, 1981; Royrvik and Smith, 
1984). Some radar investigations have indicated very narrow angular spectra 
for layers below 75 km (Fukao et al., 1980), suggesting very flat structures 
with horizontal scales greater than one Fresnel zone (Hocking and Vincent, 
1982), making turbulence an unlikely mechanism for generating such layers. 

In the absence of coincident rocket and radar observations and under the 
limitations of the one-dimensional nature of rocket measuranents , rocket data 
cannot discern the exact three-dimensional nature of the layers. One can 
only determine the consistency of rocket data with theoretical predictions. 
This chapter discusses theories f homogeneous, isotropic turbulence, 
providing a base for comparison with experimental results to be presented in 
Chapter 5. A complete development of turbulence is not intended in this 
chapter. More extensive presentations are available by Tennekes and Lumley 
(1972) and Lumley and Panofsky (1964). 

Section 3.2 gives a qualitative description of turbulence and its 
associated parameters ar.d measurable features. Section 3*3 reviews various 
sources of atmospheric turbulence. Section 3.4 discusses the Navier-Stokes 
equations, the problem of their closure, and their spectral interpretation. 
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Dimensional and similarity arguments are presented in Section 3.5, leading 
to a power law spectrum for certain length scales. Various theories for the 
spectrum at high wave numbers are examined in Section 3.6. Also, a scalar 
spectral model by Hill and Bowhill (1976) is applied to predict the type of 
spectra expected from the rocket data in this thesis. 

3.2 Qualitative Description of Turbulence 

Turbulence is often described as a cascade of energy from larger scales 
to smaller scales. A turbulent flow may be visualized as consisting of 
eddies of different sizes. The various scales of whirlpools and eddies seen 
in a river downstream from a bridge pier provide a graphic example of 
turbulence. These scales may be divided into three regimes with respect to 
their mechanisms of energy transfer. The first regime contains the largest 
scales (e.g., size of the bridge pier) in the flow. They draw their energy 
from instabilities in the mean flow. All smaller eddies find themselves in 
the rate of strain field of these large eddies, and energy is transferred to 
the smaller scales through vortex stretching (Tennekea and Lumley, 1972). 
Tennekes and Lumley also show that most of the energy transferred to an eddy 
comes from eddies with only slightly larger scales. Thus an energy cascade 
develops, energy being continuously transferred to decreasing scale sizes. 

In the second regime the mean flow no longer contributes energy to the 
turbulent flow. Bnergy is transferred to the second regime only through 
cascading from the larger eddies in the first regime. The cascading process 
continues, with little energy lost through ether processes. This regime is 
often called the inertial regime. The cascading energy eventually is 
transferred to the eddies of the third regime, containing the smallest 
scales in the turbulent flow. At these small scales viscosity removes 
significant amounts of energy from the cascade. Decreasingly smaller scales 
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receive less and less energy ae it is ranoved by viscous dissipation. This 
regime is often called the viscous range. 

Fluids and fluid flows may be described by a variety of parameters. 
Dynamic viscosity p, a constant for a fluid at a given temperature and 
independent of pressure, gives the proportionality between a fluid flow's 
shear stress and rate of shear strain. Kinematic viscosity, v, a more 
common term in describing fluids, is the dynamic viscosity divided by the 
density of the fluid. Due to this dependence on density, the kinematic 
viscosity is a strong function of pressure, in contrast to the dynamic 
viscosity. Viscosity damps velocity shears, dissipating kinetic energy as 
heat. D, the molecular diffusiyity, describes the rate of diffusion of a 
scalar contaminant in the fluid. 

~ + V-VN = DV 2 N (3.1) 

3 1 ~ - 


where N is the scalar concentration and V is the velocity field. In analogy 
to v, a fluid with a high D rapidly dissipates any large gradient in the 
scalar concentration. Another useful parameter, the Schmidt number, may be 
defined from v and D: 


D 

Substitution of the thermal diffusivity, <, for D gives the Prandtl number. 

A fluid flow, described by a velocity field, may be characterized by 
another set of parameters. The Reynolds number ie given by 


Re = 


VL 


(3.2) 


where V is a characteristic velocity of the flow, L is a characteristic 
length, and v is the kinematic viscosity of the fluid. For flows with 
Reynolds numbers higher than some critical value, the flow becomes unstable. 
This number must be found experimentally and depends on the boundaries of 
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the flow. For stratified flows, the Richardson number, Ri, provides another 

criterion for the creation of turbulence: 

2 

Ri - - r (3.3a) 

OV /3z) 
o 

where 

a, 2 = -(L-l)g 2 + (3 ,3b) 

B 2 t 3z 

c 

is the Brunt-Vaieala frequency, V Q is the horizontal flow's mean 
velocity, z is perpendicular to the planes of stratified flow, X is the 
ratio of the specific heats, g is the acceleration due to gravity, c is the 
speed of sound, and T is the temperature. The negative of the temperature 
gradient, - is known as the lapse rate. Note that for sufficiently 
high lapse rates the Brunt- Vaisa la frequency becomes imaginary. The lapse 
rate at which w « 0 is known as the adiabatic lapse rate. An atmosphere 

D 

in adiabatic equilibrium assumes this lapse rate. A lapse rate greater than 
the adiabatic lapse rate is termed a superadiabatic lapse rate. If Ri is 
less than unity the decay rate of turbulence is smaller than the generation 
rate. Under such conditions a turbulent flow is sustained (Hodges, 1967) . 
However, Ri less than unity is not a sufficient condition for the creation 
of turbulence in a laminar flow. For this to occur Ri must fall below some 
critical value Ri . Using theoretical arguments Miles (1961) and Howard 
(1961) propose a value of 1/4. 

3.3 Sources of Turbulence in the Atmosphere 

A number of mechanisms are believed responsible for the creation of 
atmospheric instabilities necessary for the generation of turbulence. Two 
important instabilities discussed here are the convective instability and 
the Kelvin-Helmholtz instability (KHI) , Convective instability is produced 
by lapse rates greater than the adiabatic lapse rate. Such an inversion can 
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be produced by large amplitude gravity waves altering the background 
atmosphere. The KHI is produced in unstable wind shears characterized by a 
Richardson number less than Ri cr . This instability may also be induced in 
a stable wind shear by interaction of a gravity wave and a critical layer 
within the shear. A detailed theoretical treatment of the two instabilities 
is available in Chandrasekhar (1961). Measurements of the instabilities in 
the atmosphere are reviewed by Gossard and Yeh (1980). 

Convective, or Rayleigh-Taylor , instabilities are the result of a 
vertical temperature gradient smaller than that of an adiabatic atmosphere. 
If a parcel of air at position z and temperature T(z) is displaced 
adiabatically and vertically by Az it expands, acquiring a temperature T(z) 
+r Az where r B -g/c_, g is the acceleration due to gravity at height 
z, and 0 ^ is the specific heat of air at constant pressure. However, the 
ambient temperature at z + Az is T(z + Az) . If this temperature is greater 
than T(z) + r Az, then the parcel is cooler and hence heavier than the 

o 

ambient air. The parcel tends to return to its initial position. Overshoot 
and subsequent oscillation of the parcel may lead to the excitation of a 
gravity wave. Conversely, if the ambient temperature is less than T(z) + 
r Az, the parcel continues to rise. The convection of these parcels 

fl 

creates turbulent eddies as the parcels move through the ambient atmosphere. 

Convective instability is more generally described as the result of an 
2 

atmosphere in which to- is negative. In such a situation m is 

B B 

imaginary, and a displaced air parcel does not oscillate but continues to 
accelerate away from its initial position. A negative to 2 also forces 
the Richardson number below zero. A negative Richardson number thus 
indicates convective instability. The Brunt-Vaisala frequency may be more 
conveniently defined in terms of potential temperature. The potential 
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temperature, 6, of an air parcel ie that temperature the parcel would assume 
if adiabatically compressed or expanded to a pressure of 1000 mb. 


0 = T(- 


1000 


1=1 

) V 


where P, the pressure, is given in millibars. The Brunt-Vaisala frequency, 

in terms of potential temperature, is 

2 _ £ j)B 


8 3z 
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The lapse rate at which - r— “ 0 is the adiabatic lapse rate. An 

d Z 


atmosphere in adiabatic equilibrium has 
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7z 


0 for all z. Under such a 


condition a displaced air parcel neither oscillates nor accelerates. If 
- t — is greater than zero, a superadiabatic lapse, then a is 

3Z n 

negative. Thus the sign of the potential temperature gradient determines 
stability. 

Such convective instabilities can result' from gravity wave growth. As 
a gravity wave propagates upward' its amplitude increases with decreasing 
density to preserve energy flux. The wave eventually reaches an altitude 
where the increased amplitude significantly modifies the background 
atmospheric parameters, producing regions of superadiabatic lapse rate. 
Hodges (1967) calculates the Richardson number as modified by a 
monochromatic gravity wave and finds negative values at periodic intervals 
for sufficiently large wave amplitudes. These planar regions of convective 
instability propagate with the wave, are normal to the direction of wave 
propagation, and occur at the minima of the wave's wind shear. 

The KHI is produced in regions of high vertical shear in a horizontally 


stratified flow. If Ri falls below the critical value, it is energetically 
feasible for vertical perturbations to overcome the stabilizing effect of 

i 

decreasing density with altitude. The perturbations grow and the vorticity 
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accumulates in billows forming a wave pattern, the billows continue to draw 
vorticity from the mean flow, eventually forming what are often referred to 
as cat's eyes. The growth of the cat's eyes from the billows is numerically 
simulated by Patniak et al . (1976). The eventual breakdown of the cat's 
eyes produces turbulence. Gossard et al, (1970), using a high resolution 
radar, show the development of KHI in the lower atmosphere and the resulting 
turbulence. 

An initially stable flow can develop KHI through interaction of a 
gravity wave and a critical layer in the flow. A critical layer exists in a 
shear at that level where the flow's velocity equals the horizontal 
component of the incident gravity wave's phase velocity (Booker and 
Bretherton, 1967). The analysis of Booker and Bretherton (1967) models the 
interaction using a monochromatic wave incident upon an inviscid, adiabatic, 
Boussinesq fluid. The flow has an initial Richardson number everywhere 
greater than 1/4. Booker and Bretherton (1967) find that the gravity wave 
is strongly attenuated upon passing through a critical layer and that the 
wave's energy is deposited into the flow without loss through the generation 
of turbulence or other dissipative processes. Geller et al . (1975), using a 
similar model, show that the deposited energy modifies the flow such that 
regions of dynamic instability are produced. However, the model is 
incapable of showing the detailed initial development of the instability. 

Such a detailed simulation is described by Fritts (1978). Neither 
simulation can demonstrate the final turbulent breakdown of the Kelvin- 
Helmholtz instabilities . 

3.4 Statistical and Spectral Approach 

So far only a qualitative description of turbulence has been given. 

Due to the chaotic and three-dimensional nature of turbulence, a mathematical 


25 




treatment requires the use of random variables and statistical methods. 
Despite this statistical treatment the Navier-Stokes equations for fluids 
have no closed form solution, requiring the addition of observational and 
physical arguments. Such argument. are most easily applied to the power 
spectra of the velocity and advected scalars. Simpler spectra may be 
derived with various statistical assumptions. This section shows the 
development of these spectra and the restrictions on their use. 

As with the kinetic theory of gases, turbulence requires a statistical 
approach. Although the motion of a fluid in a turbulent velocity field is 
in theory a classically deterministic process, the amount of information 
required to treat it as such makes a deterministic approach wholly 
impractical. The statistical treatment essentially examines a quantity's 
fluctuation about a mean and the spatial dependence of averages of those 
fluctuations. The quantity might be a velocity component, a temperature, or 
the density of an advected scalar in a velocity field. For example, the i'th 
velocity component V^(x,t) and a scalar's density N(x»t) can be decomposed 
into a mean and a fluctuation: 

V^Xjt) ■ ^( 5,0 + u^ (x , t ) (3.4a) 

N(x,t) * N(x,t ) + n(x,t) (3.4b) 

where 

D i (x,t) ■ <V.(x,t)> 

N(x,t) » <N(x,t)> 

and < > denotes an ensemble average. The calculated averages of most 
interest are the covariances of the fluctuations. These may be written as 
<u i (x,t)u j (x , ,t)> 

<n(x,t)n(x’ ,t)> 


(3.5a) 

(3.5b) 




r<V*» * -W-. j. 
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To make practical measuranents of the covariances bob’? assumptions must 
be made. First, a moving probe (e.g., rocket borne Langmuir probe) is often 
used to measure scalar density along a path. Due to the finite velocity of 
the probe the turbulent field evolves during the time of the measuranent. 

To make meaningful spatial covariances the turbulence must be assumed to be 
"frozen" during this time. For example, a rocket travelling at 1 km/sec 
traverses a turbulent layer 1 km thick in 1 second. Radar measurements 
(Rottger et al., 1979) indicate that turbulent layers have lifetimes on the 
order of minutes. The frozen-in assumption is hence quite applicable to 
rocket measurements of such layers. A second assumption, ergodicity, must 
be made from a practical point of view. This allows for replacement of the 
ensemble averages by spatial or time averages. In the rocket experiment 
described above, spatial averages are used. As a different example, 
measuranents of wind turbulence using a fixed, ground based hot wire 
anemometer would replace the ensemble average by a time average. 

Under the conditions of statistical homogeneity and statistical 
isotropy the covariances may be simplified. Statistical homogeneity, 
hereafter referred to as just homogeneity, implies that the probability 
density of the fluctuations is not a function of position. In addition, the 
averages <V^> and <N> do not depend on position. Under this condition 
<U£(x,t)uj(x't)> can only be a function of r ° x' * x. Due to the 
inherent inhomogeneity of the large energy producing scales, such an 
assumption can only approximate an actual flow (Batchelor, 1953), and the 
less restrictive, more realistic, assumption of local homogeneity is 
necessary. Under this condition homogeneity is restricted to small-scale 
variations. The covariance <u. (x,t)u •(x'’ ,t)> is then independent of 
position only for small r. Statistical isotropy, a refinement of 
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statistical homogeneity (Kolmogorov, 1941), implies that there is no 
preferred direction in the velocity or scalar field. As with homogeneity, 
the qualifier "statistical 11 is implied. It is quite conceivable that a 
field be homogeneous yet anisotropic. For example, a field may have a 
homogeneous preference for u^ components over Uj components, thus 
representing a homogeneous anisotropic field. Isotropy, like homogeneity, 
is most realistically applied at small scales. This condition is called 
local isotropy. When both isotropy and homogeneity hold, the covariances of 
and n depend only on [x - x'|. Kolmogorov (1941) claims that for flows 
with high Reynolds number a wide gap in k-space exists between the large and 
small eddies. Under such conditions Kolmogorov hypothesizes local isotropy 
for the smaller scales. Tennekes and Lumley (1972) make a similar claim by 
examining the interdependence of different sized eddies. They find that 
although the strain rate field of a large eddy induces anisotropy in passing 
energy to a small eddy, at sufficiently small Beales the anisotropy is 
temporary, and the small eddy's mean state is isotropic. 

Fluid flows may be described by the Navier-Stokes equations; 

8V 12 /v 

— + (w)V = - — VP + vV V (3.6a) 

9 1 - - - p - - ~ 

vv - o (3 - 6b) 

•V M 

where V is the velocity field, P is the pressure, and p i® the fluid 
density. Equation (3.6b) states that the fluid is incompressible. As 
mentioned earlier, turbulent flows must be described statistically. 
Decomposing V. and P into U. + u. and P + p evaluated at x and t, 
assuming the tensor summation notation, and applying the incompressibility 
condition gives 


9(U i +U i ) 3(U t +u I )<U dt +u,) 


9t 


* i i r = i icg±E > , 


9x 




9x. 


a (u i +u i ) 

9x. 9x. 

3 3 
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Multiplying by u„(x',t) and averaging gives an equation in terms of 
9 - 

<u„(x',t)u,(x,t)>, <u„(x't)u 4 (x,t)u.(x,t)>, and other terms. 

Attempting to solve for the triple correlation using the Navier-Stokes 
equations leads to an equation involving a fourfold correlation. ThiB 
continued predicament is referred to as the closure problem. Various 
statistical assumptions may be made to allow for an approximate numeric/il 
solution. Some of these methods are reviewed by Tatsumi et al. (1978). 

Hill and Bowhill (1976) point out that even under isotropy and homogeneity 
the equations only reduce to a single equation involving two unknowns: the 
double and triple correlations. It should also be noted that the diffusion 
equation 

f? + v*vn = dv 2 n 


is coupled to the Navier-Stokes equations, resulting in a closure problem 
for scalar turbulence also. 

It is often physically enlightening to take the spatial Fourier 
transform of the covariances. The resulting spectrum facilitates comparison 
of the energy at various scales. The following spectral formulae are taken 
from Hill and Bowhill (1976). The three-dimensional velocity fluctuation 
power spectrum, assuming homogeneity, is given by the Fourier transform of 
the covariance: 


(2it) 


-j jjj <u(x,t)»u(x + r.t)> exp(-ik'r) 


or 


Physically, this gives the strength of a velocity fluctuation at the wave 
vector k. Since | k | = 2 it/ X the spectrum also is a measure of fluctuations 
with scale size X. This spectrum can be averaged over all directions to 
give the three-dimensional energy spectrum, commonly referred to as just the 
energy spectrum: 
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E(k, t) 


i \\ k2 


— <u(x,t)*u(x + r,t)> exp(-ik*r) dr 
<2ir) J 

(3.7) 


where df^ ib the solid angle in k space and k ■ |k|. The factor of ,1/2 
normalizes the spectrum such that E(k,t) integrated over all k gives (1/2) 
<u*u>, called the turbulent kinetic energy. For a given magnitude k, 
E(k,t) gives the average energy per unit mass per unit wave number. An 
analogous three-dimensional spectrum also may be constructed for scalar 
quantities. Assuming homogeneity, the three-dimensional scalar power 
spectrum and its spatial average, the scalar spectrum, can be written as 


S(k, t) 

r(k,t) 


HI 


dr 


(2ir)‘ 


<n(x, t)n(x+r,t)>exp(-ik*r) (3.8a) 


k 2 drt, 


(2tt) 


^ <n(x,t)n(x + r,t)> exp (-Ik* r) dr 

(3.8b) 


Like the velocity energy spectrum, T(k,t) is an average over all directions 
of k. r(k,t) is normalized so that 


<n 2 > 


r(k,t) dk 


0 


Physically, f(k,t) gives the strength of the scalar magnitude's fluctuation 
at a given wave number, whereas S(k,t) gives the same for a particular wave 
vector. The units of P(k,t) are mean squared density per unit wave number. 
The two types of spectra, three-dimensional scalar power spectrum and its 
spatial average, the scalar spectrum, serve different functions in studying 
turbulence. The three-dimensional scalar power spectrum plays a major role 
in the theory of coherent scatter of radiation from atmospheric turbulence. 
The scalar spectrum is for the most, part a theoretical tool for describing 
spectral energy in a field. Yet a third spectrum is needed to interpret in 
t-itu experimental data. 
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Many in situ measurements of turbulence are made in only one 
dimension with a rapidly moving probe. The turbulent field's spatial 
variation may be derived from the probe's temporal signal and speed relative 
to the field. However, for a measured variation of k^ , wave vectors k 
with |k| > kj and k not parallel to the probe's path alias as k^ 
variations. One-dimensional spectra are defined as follows: 


or 


“Kk^, t) 
4> Ck^, t) 


dx <u (x,y,z,t)u (x+x ,y,z,t)>exp(-ik 1 x ) 
IX XX XX 
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,y,z,t)>exp(-ik 1 x 1 ) 
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i f 00 

- j d* 1 n<x+x 1 ,y,z,t) expC-ik^) 


(3.9b) 


for the scalar spectrum. The second forms of equations (3.9a) and (3.9b) 
are usually chosen for computational ease. In the isotropic case these are 

related to the energy spectra by 

>00 

**!.*> - j 20^1 (1 . k 2 /k 2 > dk O.lOa) 

k l 

,oo 

VOtj^.O » J ~-- k - » t) dk (3.10b) 

k l 

^ is the k space sum of r(k,t) over all k > kp To obtain VCk^t) each 
r(k,t) value is weighted by the ratio of the wave vectors from k to k + dk 
contributing to the aliasing and the total number of wave vectors ending on 
a sphere of radius k. The energy relationship, equation (3.10a), is given 
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by Bill and Bowhill (1976). The above equations are valid only under 
isotropy. Consequently if the turbulent field is locally isotropic only for 
kj > k|', then the equations apply only to that range. It should be 
noted that for r(k,t) or E(k,t) proportional to a power of k, the related 
one-dimensional spectrum has the same charateristic power law. If such a 
power law is limited to a particular range of k, then the observed break 
points in the one-dimensional Bpectrum differ from those in the energy 
spectrum. Figure 3.1 compares a scalar Bpectrum and its corresponding 
one-dimensional spectrum calculated from equation (3.10b). 

Under the restrictions of homogeneity and isotropy , the Navier-Stokes 
equations yield the following for E(k,t) and T(k,t) (Hinze 1975): 

8E(kj _ t 2 B T (k t ) = _2vk 2 E(k,t) (3.11) 

at u 

9r - T(k,t) = -2Dk 2 r(k,t) (3.12) 

o C , 

where T u (k) and T(k) are called the spectral transfer functions. T u (k) 
has the units energy per unit mass per unit wave number per unit time and 
represents the loss of energy by viscous dissipation at c. scale size k. 
Metaphorically, it is the leakage in the energy cascade described in 
Section 3.2. T(k) has the units of mean squared scalar fluctuation per unit 
wave number per unit time. It too represents dissipation of spectral 
content at a wave number of size k. Two important parameters may be defined 
from the spectra r(k,t) and E(k,t). These are the energy dissipation rate 
and the mean squared scalar dissipation rate x ; 

X - 2D f k 2 r(k,t)dk 

e • 2v j k 2 E(k,t)dk 
'0 

e measures the rate at which kinetic energy is dissipated by viscosity. 




Similarly, x measures the rate at which mean squared fluctuations are 
dissipated by molecular diffusion. The.-’e quantities aid in the dimensional 
analysis of the turbulence spectrum. 

3.5 Dimensional Analysis and Turbulence Spectra 

Applying local homogeneity and isotropy, dimensional and similarity 
arguments leads to a power law form of the energy and scalar spectra in the 
inertial range. The range over which this dependence applies is derived in 
a similar fashion. Spectra derived by Batchelor (1959) and Batchelor et al. 
(1959) for wave numbers above the inertial range are also presented here. 
Models of a scalar flux function, defined from the spectral transfer 
function T(k), are examined and lead to refinements of the simple power law 
spectra. 

In addition to local isotropy and homogeneity Kolmogorov (1941) 
introduced r third concept: local similarity. Under local similarity, if 
them is a large wave number range between the scale sizes at which the 
curbulence is created and the scale sizes at which viscosity becomes 
important, then the spectrum for large k is independent of the detailed 
morphology of the large scales at which energy is drawn from the mean flow. 
Such a wide range of scales is present in flows with high Reynolds number 
(Kolmogorov 1941). Due to this independence turbulence produced by 
different sources in different media is found to be comparable at scales 
where local similarity holds. The comparison is subject only to certain 
parameters of the fluid and the flow. For example, a knowledge of these 
similarity parameters allows for comparison between turbulence generated in 
the ocean and turbulence generated in the atmosphere although the media, 
energy sources, and scales are quite different. Local similarity is obeyed 
under local statistical equilibrium (Batchelor, 1953). This requires that 


the time evolution at email scales proceed much more rapidly than at large 
scales. Equilibrium is established at small scales with negligible influence 
by the large-scale evolution. Small-scale structures may then be studied in 
general without considering the detailed time evolution at the large scales 
or the particular method of the turbulence's generation. 

Before deriving the various spectra, a few characteristic wave numbers 
need to be defined. As mentioned previously, turbulent energy cascades from 
larger scales to smaller scales, the scales being divisible into three 
regimes with respect to their methods of energy transfer. The concept of an 
inertial range may be extended to scalar fluctuations. The analogous regime 
is called the convective range: Bcalar fluctuations therein are not 
dissipated by diffusion. Similarly, the diffusive range is analogous to the 
viscous range, k^, called the inner scale or Kolmogorov scale, is a 

measure of the wave number at which viscosity begins to significantly damp 

the velocity fluctuations, identifies the boundary between the 
inertial and viscous ranges. k b , sometimes called the Batchelor wave 
number, defines a scale beyond which diffusion damps out the scalar 
concentration fluctuations which are already subject to viscous damping, 
k , referred to as the Obukhov-Corrsin wave number, defines a scale beyond 
which diffusion effects are important in the absence of viscous damping. 

k d - (s/v 3 ) 1/4 (3.13a) 

k b » ( e /vD 2 ) 1/4 (3.13b) 

k c ■» (e/D 3 ) 1/4 (3.13c) 

Batchelor (1959) shows that k c applies only when D » v and that kj, 
applies when v » D. It should be noted that scalar fluctuations may 
persist for wave numbers beyond k d when v » d. In that case k b » k^. 

Assuming similarity, isotropy, and homogeneity to apply locally. 


Kolmogorov (1941) hypothesized that the structure of turbulence should 
depend only on the similarity parameters e and v. Using dimensional 
analysis Tennekes and Lumley (1972) find that 
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E(k) - e m v n f(k/k d ) 

where m and n are constants, f(k/k d ) is a nondimens ional function, and 
k^ depends only on e and v. Given the units for v, e, and E(k) the only 
possible choice for m and n is 1/4 and 5/4, respectively . This equation 
implies that for all flows E(k)/(ev^)^^ should yield the same 
dimensionless curve f(k/k d ) at scales subject to homogeneity, isotropy, 
and similarity. Kolmogorov (1941) further argued that in the inertial range 
viscosity is unimportant. E(k) then depends only on e. This condition 
requires that f(k/k d ) cancel out the constant in front, leading to 

the energy spectrum 

E(k) ■ k « k d 

where a is a nondimensional constant called the Kolmogorov constant, k « 
k d is a reminder that this spectrum applies only for k for which viscosity 
is unimportant. It should also be noted that the above development assumes 
locally a steady state turbulence field, eliminating the t dependence from 
E(k,t). 

By similar methods Boston and Burling (1972) use the similarity 
parameters X, v, e, and D to find the form for the scalar spectrum: 

r(k) = xe" 3/i!t v 5/ S(v/D,k/k d ) 

where t«V/D,k/k.) is a dimensionless function. An inertial range may be 

found by disallowing effects of v and D. 

r(k) = gye ^/\”^/3 k << smaller of k, and k 

d c 

where 0 is a dimensionless constant. 

The shape of the scalar spectrum beyond the Inertial range has been 
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studied by Batchelor (1959) and Batchelor et al. (1959). The relation 

between D and v affects not only the cutoff wave numbers k^, and 

k but also the spectral shape and slope. For D « v the scalar 
c 

fluctuations can occur at wave numbers beyond k^ in spite of the loss of 
kinetic energy through viscous damping. Batchelor (1959) shows that this 
leads to a scalar spectrum 

r(k) = - xo 1 exp(Dk‘'/cj) (3.14) 

where a is an average strain rate estimated by Batchelor (1959) to be 

a = - .5(e/v) 1/2 (3.15) 

The argument of the exponent in equation (3.14) then becomes -2(k/k^). 

Thus for k d « k « k fa the scalar spectrum varies as k This range 
is termed the viscous-convective range since viscosity dominates the energy 
spectrum and diffusion is unimportant in the scalar spectrum. The range k > 
k. is called the viscous-diffusive range since both viscosity and diffusion 

D 

affect the spectra. In the case D » v Batchelor et al. (1959) predict 

r(k) = (1/3) a X e 2//3 D k « k << k,. 

c a 

This range is called the inertial-diffusive range since the scalar 
fluctuations are damped but viscosity has little effect. The upper k range 
where both viscosity and diffusion become important is called the 
viscous-diffusive range just as in the D « v case. These ranges are 
illustrated in Figure 3.2. The middle spectrum gives a hypothetical energy 
spectrum. The expected scalar spectra for v » D and v « D are shown above 
and below the energy spectrum. For D » v the viscous-diffusive range is 
not shown. This range is not considered by Batchelor et al. (1959). Due to 
the very steep slope in the inertial-diffusive range extremely low values of 
T(k) i f/re expected in the viscous-diffusive range. Relevant measurements and 
studies are not abundant. In the case D % v the viscous and diffusive 
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effects occur at approximately the same wave number. This situation is not 
addressed by Batchelor (1959) or Batchelor et al. (1959). Its discussion 
here is deferred until the presentation of the model by Hill and Bowhill 
(1976) in Section 3.6. 

As an alternative to directly estimating r(k) , the spectra may be 
derived by modelling T(k) and working from the Bteady state equation (3.12). 
Models of T(k) leading to scalar spectra are extensively reviewed in Hill 
and Bowhill (1976). Three of the models are presented here. To model T(k), 
the spectral flux function F(k) needs to be developed. 

- T(k) (3.16) 

3k 

Physically, F(k) describes the rate of the spectral energy cascade through a 
scale size with wave number k. Its units are mean squared scalar per unit 
time. In the inertial range the flux should be approximately constant and 
equal to X. Since viscosity is unimportant} T(k) should be quite small. 

T(k) should increase considerably in the diffusive regime. 

Using results for spectral flux from Corrsin (1964), Pao (1964) shows 

that 

T(k) - k" 5/3 exp[(-3/22,)(k/k c ) 4/3 ] k « kj (3.17a) 

r(k) « k 1 exp[ (-D/&'o)k 2 ] k >> k^ (3.17b) 

where & and if are nondimensional constants and a is the strain rate given 
by equation (3.15). If v » D and a is assumed to be of the order of unity, 
then k fi » k d , making the exponential term in equation (3.17a) close to 
unity for k « k^. This results in the familiar -5/3 power law for the 
inertial range. For -1 equation (3.17b) is of the Bame form as Batch- 

elor's, equation (3.14), predicting an exponential decrease with k above 
k b . If, however, D » v then equation (3.17a) predicts a rapid decrease in 
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in r(k) for k « k c * This decrease with k is steeper than the -17/3 power 
law predicted by Batchelor. 

For high Pr or Sc in the viscous-convective and viscous-diffusive 
ranges a model by Kraichnan (1968) predicts the scalar flux function to be 

F (k) - ± k 4 fk" 2 r(k)] 

where A is independent of wave number but increases with Reynolds number. 
Assuming a steady state in equation (3.12) gives (Mjolsness, 1975) 

T(k) - 5(x/A)k _1 (l+Z)exp(-Z) 

where 

Z = (30D/A) 1/2 k 

Mjolsness (1975) also derives the corresponding one-dimensional spectrum 
V(k 1 ) = 5( X /A) k^ 1 exp(-Z^) 

If Z « 1 the scalar spectrum is approximately 
r(k) ■ k _1 

the dependence found by Batchelor (1959) in equation (3.14). For large k, Z 
also becomes large, causing an exponential decrease in the spectrum. 

Hill and Bowhill (1976) also analyze a model by Howells (1960) which is 
applicable for all Pr and Sc. They find that the proportionality 

r(k) « k" 5/3 [c + c * ( — ) -8/3 i 

k 1 

c 

holds in the inertial- convective range with C and C' constants. C' is taken 
to be small since the spectrum should approach a -5/3 power law for k « 

^ as k becomes small although still in the inertial range. The viscous- 
convective and viscous-diffusive range spectra agree with the form of 
Batchelor 's, equation (3.14), For D » v the model predicts an inertial- 
diffusive spectrum of the same form as equation (3.17a). 


In summary, all models agree on the -5/3 power law for k well inside 
the inertial range. For v » D most models predict a scalar spectrum pro- 
portional to k -1 for k d « k « k b . The models differ substantially 
for the inertial-diffusive range (D « v), Batchelor et ai. (1959) 
predicting a power low and other models predicting an exponential decrease 
above k fi . None make solid claims to the spectral behavior for D ^ v. 

This is unfortunate as such conditions are found in the middle atmosphere. 

3 .6 Estimation of Spectral Form in the Middle Atmosphere Using a Model by 
Hill and Bowhill (1976) 

The scalar spectrum models given so far say little about the spectrum 
between ranges and apply only for the asymptotic cases of v » D and v « D. 
Hill and Bowhill (1976) have devised a spectral model, based on Corrsin 
(1964) and Pao's (1964) results for F(k) to describe the spectrum for all 
spectral ranges and all values of Sc and Pr. This section presents Hill and 
Bowhill's (1976) model with some clarifications added by Hill (1978). A 
spectral shape for middle atmospheric turbulence, based on estimates of Sc 
for those altitudes, is also predicted. 

The Corrsin-Pao model determines the scalar spectrum by modelling the 
scalar flux as (Corrsin, 1964): 

F(k) * s(k) (k) 

where s(k), called the flux rate, is estimated by Pao (1964) to be 

s(k) - Z e 1/3 k 5/3 k « k d (3.18a) 

s(k) » Z'ok k » k d . (3.18b) 

This applies to both the cases D » v and v >> D (Pao, 1964). Pao describes 
s(k) physically as the rate that a spectral element is tranf erred across k. 
Since r (k) is a measure of spectral content at k the product of T (k) and 
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e(k) is dimensionally the spectral flux F(k), as Corrsin'B hypothesis 
states. Hill and Bovhill (1976) hypothesise that the flux rate dependence 
on k varies from 5/3 to 1 in the k range close to kj. The scalar spectrum 
then also has a smooth transition. A characteristic wave number k* located 
between the inertial-convective and viscous-convective ranges may be 
defined; 


3 lnfr(k)] „ _ A 

3 tln(k) J k 0 k * 3 


for V >> D 


Although the corresponding one-dimensional scalar spectrum has a shape simi- 
lar to the three-dimensional scalar spectrum, it is found that kj* < k* 
where k^* is the transitional wave number for the corresponding one-dimen- 
sional spectrum. To form s(k), and subsequently F(k), Hill and Bowhill 
(1976) form a nondimensional wave number expression, 


K = ln(k/k*) 

(3.19) 

and create a parameter 

4> (K) = -4/3 +1/3 tanh(aK) 

(3.20) 


which varies from -5/3 to -1 over the range of all k. The value of a deter- 
mines the rate at which <KK) changes from -5/3 to -1, s(k) is then proposed 
to be 

s(k) = s(k Q ) ex P 

where kg is the lowest wave number in the inertial-convective range. To 
find r(k), F(k) is formed by multiplying s(k) and T(k) and substituting into 
the steady state version of equation (3.12), giving a differential equation 

for r(k). Taking kg to approach zero. Hill and Bowhill (1976) find the 
solution to be 



r(k) - X0E -1/3 (k *)“ 5/3 r(K) 
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where r(K) ie given by 


ln[r(K)J = - ^ K + ~ ln[2 cosh(aK)] 


- C 


K e 5y/3 [2 cosh(ay)] l/3a 


dy 


£ » 20(k*/k ) 4/3 
c 

k*/k. and a are assumed constant over all v/D (Hill, 1978). The model 
a 

makes no assumptions about D's relation to v. Although k* is defined above 
for D « v, Hill (1978) points out that for D » k* represents the 
transition wave number between the inertial-diffusive and the viscous-diffu- 
sive range. Figure 3.3 shows the calculated normalized one-dimensional 
spectra as a function of k^/k* for various values of the parameter £ where 
symbols with subscript 1 are associated with the corresponding one-dimen- 
sional spectrum fOt^). 

It is of interest to estimate the type of spectra the above model 
predicts for middle atmospheric parameters. Since only a general shape is 
desired, the important parameter to estimate is £. Many experiments have 
made measurements of the constant Pj. . Experiments reviewed by Hill and 
Bowhill (1976) give a range centered about 0.6, They estimate that 3 « 

5^j/3, making 3 of the order of unity. Values of k c and k d vary with 
altitude. Three parameters are necessary to compute k £ and k^: e, 

and D, As an example, at. 80 km D is 1 .0 x 10 cm /sec (Heicklen, 1976), 
and v is 7.15 x 10 m /sec (U.S. Standard Atmosphere, 1976). Since v is 
so close to D (i.e., Sc % 1), k^ and k c are quite close together, k*, 
the transition wave number between the inertial- and viscous-diffusive ranges 
is between these two numbers, making k*/k c close to unity. Thus £ is of 
the order of unity, independent of e. Referring to Figure 3.3 it is observed 
that Hill's model predicts neither a well-defined inertial-diffusive or 




Figure 3.3 One-dimensional spectra for various values of the parameter £ (Hill and Bovhill, 1976). 
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viecouB-convective range. The elope of the spectrum should be -S/3 up until 
a rolloff at the characteristic k value of k d or k c , whichever is smaller. 
Comparisons to actual data from the University of Illinois rockets are given 
in Chapter 5. 


4. FINE STRUCTURE EXPERIMENT AND METHODS OF DATA ANALYSIS 


4.1 Introduction 

To measure the spectrum of electron density irregularities in the D 

region a Langmuir probe is mounted on a sounding rocket and flown through 

the ionosphere. The fine structure experiment flown by the University of 

Illinois is described by Klaus and Smith (1978). All flights have carried 

this same configuration, with two exceptions described in Roth (1982b). In 

the fine structure experiment the probe current, assumed proportional to the 

electron density over regions of constant electron temperature, is measured 

as a function of time. Rapid current fluctuations superposed upon the 

relatively slowly varying mean current are sometimes observed. Signal 

processing on board the rocket converts £;hese current fluctuations to 

AN /<N > where AN is the fluctuation in electron density about a mean 
e e e 

Value <N g >. This signal is then telanetered to ground and recorded on an 
FM tape recorder. Analog preprocessing of these data helps to distinguish 
possible ionospheric irregularities from interference caused by other 
experiments on board the rocket. 

Once suitable regions of interest have been established, the analog data 
are digitized for further processing. The digitizer consists of an 8-bit A/D 
peripheral card controlled by an Apple II microcomputer. The data are 
digitized at a 5 kHz rate and stored on floppy disk. Each data file con- 
tains 25000 bytes, or 5 seconds of data. These data are then transferred 
to a CDC Cyber mainframe computer. After additional formatting the data are 
segmented and Fourier analyzed with a FORTRAN program called SPECTRA. This 
program estimates the power spectrum of a data segment by calculating the 
square magnitude of the discrete Fourier transform (DFT) of the segment. 
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It should be emphasized that these spectra represent a one-dimensional 
spatial spectrum, where isotropy in assumed in converting the temporal 
variations to spatial variations. The theories of turbulence described in 
Chapter 3 are concerned with three-dimensional spectra, and theoretical 
one-dimensional spectra are derived therefrom only under the assumption of 
isotropy. Since no independent measure of isotropy iB usually made during 
the rocket experiments, spectra calculated in this chapter cannot be 
directly compared with the theoretical one-dimensional spectra; they can at 
most demonstrate consistency with theory. 

Section 4.2 describes the University of Illinois fine structure 
experiment. Section 4.3 discusses three possible sources of noise and 
distortion. An analog preprocessing method iB presented in Section 4.4. 
Section 4.5 discusses the digitizing of the analog data. Section 4.6 
describes postprocessing on a mainframe computer. 

4.2 Fine Structure Experiment 

The fine structure experiment carried out by the University of Illinois 

uses essentially the same instrumentation as described by Klaus and Smith 

(1978). A Langmuir probe is extended from the rocket body into the 

ionosphere and held at 4.05 V above the potential of the rocket body. The 

probe may be either a single electrode at the rocket noBe or a pair of 

electrodes in parallel extended from the sides of the rocket. As the rocket 

passes through the ionosphere, spatial fluctuations in electron. 

concentration are observed as current fluctuations in time. The mean 

—8 -5 

current is usually in the range of 10 to 10 A in the D region. A 
plot of current versus height for the 28 Hay 1975 Peru flight is shown in 
Figure 4.1. The current is processed by a logarithmic electrometer and an 
ac amplifier to give the percentage fluctuation in electron density. 



Figure 4.1 Probe current from 28 May 1975 Peru flight. 
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vA, 


.. . 


AN e /<N e >, as a function of tiute. This signal, along with data from 
other experiments on the rocket, is telemetered to the ground and stored on 
an FM tape recorder. 

A block diagram of the system is shown in Figure 4.2. For fine 

structure measurements the sweep generator ia held constant at 4.05 V, the 

sweep mode being activated to collect electron temperature data. The log 

electrometer produces a voltage proportional to the logarithm of the input 

current. Its output is 0.80 V per decade of probe current. Klaus and Smith 

(1978) investigate the frequency response of the logarithmic electrometer and 

find that it acts as a low pass filter. The upper 3 dB point is determined 

by measuring the rise time of the output in response to a current input 

consisting of a square wave superposed on a dc current. At an input dc 
—8 

current of 10 A and a square wave adjusted' *ko produce a 50 mV output 
swing, Klaus and Smith (1978) find the upper 3 dB frequency to be 2.5 kHz. 
Furthermore, Zimmerman and Smith (1980) find that the 3 dB frequency 
increases with an increasing dc current component. Since currents greater 
than 10 A are expected ivr the D region (Figure 4.1), the frequency 
response of the electrometer should extend at least up to 2.5 kHz. The 
output of the electrometer is fed to an ac amplifier. This has a measured 
upper 3 dB frequency of at least 1 kHz and a gain of 100. Frequency 
response measurments vary between devices. Measurements of three amplifiers 
are presented in Appendix I. In addition, the amplifier has a single pole 
roll-off below 50 Hz, ensuring that the amplifier is not saturated by low 
frequency noise expected at the spin frequency of the rocket (4 to 8 Hz). 

To predict the output or this system it is assumed that the input 
current, I, may be decomposed into a mean and a fluctuating component: 


I - <I> + A I 



Figure 4.2 Block diagram of dc probe and fine structure experiment. 

Sweep generator is held at 4.05 V for electron density 
fine structure neasurenents (Klaus and Smith, 1978). 
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where < > represents a time average. This current is assumed proportional 
to the electron density as described in Chapter 2. The electrometer output, 

V LE’ ie 

V LE » 0.80 log ( <I> + AI ) . 

Conversion to the natural logarithm allows expansion of V LE to 


'LE 


0 80 1PA<1>) + o 80 In ( 1 AI/<I> ) 
0,80 In 10 +U,8 ° In 10 


If the fluctuations in current are small compared to the mean current, then 
the output may be approximated by 


'LE 


0.80 


In (<I» 0.80 A I_ 

In 10 In 10 <I> 


(4.1) 


The constant component of V LE disappears due to the ac coupling of the ac 
amplifier following the logarithmic electrometer: 



100 


JL80. 
la 10 


AI 

<I> 


(4.2) 


where V^ c iB the ac amplifier output. Since current and electron density 
are assured to be proportional, this is equivalent to 35 AN g /<N e >, a 
0.35 V deviation for a 1.0 % deviation in electron density. 

The data signal on the rocket is telemetered to the ground and stored 
on analog tape. Following the ac amplifier the AN /<N > signal is 

c 6 

centered about zero volts. This is stepped up to 2.5 V to conform to the 0 
to 5 V input swing allowed by the voltage controlled oscillator (VC0) of the 
telemetry system. In the IRIG FM telemetry system employed by the 
University of Illinois the VCO output frequency varies by +. 7.5 % for a 
+ 2.j V input. The modulation index is 5. Each experiment on the rocket 
is assigned an IRIG channel. As an example, Table 4.1 gives the IRIG 
channel assignment for each experiment on the 28 May 1975 Peru flight. The 
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N 


Table 4.1 IRIG telemetry channel assignments for 28 May 1975 Peru flight. 


Channel 

Number 

Center 

Frequency* 

(kHz) 

Information 

Bandwidth** 

(Hz) 

Signal 

19 

93 

1395 

RF probe, output 

18 

70 

1050 

Tip probe, fine structure 

17 

52.5 

790 

Receiver #1, modulation 

16 

40 

600 

Receiver #2, modulation 

15 

30 

450 

Tip probe, log output 

14 

22 

330 

Boom probe, log output 

13 

14.5 

220 

Boom probe, linear output 

12 

10.5 

160 

RF probe., monitor 

11 

7.35 

110 

Receiver #1, AGC 

10 

5.4 

81 

Receiver #2, AGC 

9 

3.9 

59 

Magnetometer 


‘deviation = ±7.5% 
“modulation index = 5 
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channel assignment depends on the desired bandwidth. Channel 18 is usually 
selected for the fine structure experiment, giving a bandwidth of 1050 Hz. 
The various channels are combined and transmitted to ground on a 240.2 MHz 
carrier. The telemetry ground station strips off the 240.2 MHz carrier and 
records the raw FM signal (containing all channels) on one track of an FM 
tape recorder. This track is usually called the “combiner" track. In 
addition to this data signal the ground station simultaneously records a 
number of other signals on different tracks. To aid in the digitization of 
the data a 100 kHz reference signal is recorded t i a separate track. 

Another track, called the "station multiplex" (labelled MUX), contains two 
time codes for absolute determination of elapsed time on the tape, a voice 
record of the flight progress, and miscellaneous data from instruments 
monitoring the telemetry signal. 

4.3 Sources of Noise and Distortion 

There are a number of noise sources which must be identified to 
properly analyze the data. A periodic modulation of the fine structure 
signal at the rocket spin frequency is often observed. Axial spin provides 
stability for the rocket. The spin frequency, nearly constant throughout a 
particular flight, is generally from 4 to 8 Hz. This is most accurately 
measured by an on board magnetometer measuring the component of the earth's 
magnetic field along the axis of the magnetometer. The strength and 
character of the effect on the fine structure signal varies from rocket to 
rocket, and even during the same flight the effect slowly changes with time. 
Figure 4.3 shows a segment of data exhibiting severe spin contamination. 

The most probable source of this noise is the passing of the Langmuir probe 
through the rarefied wake of the rocket. This is borne out in the 
observation that spin noise shows up strongly in flights using boom probes 
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as opposed to nose tip probes. Figure 4.3 exhibits a spin frequency of 
3.984 Hz from a rocket equipped with boom probes. The existence of two 
oscillations every spin period is consistent with the suggestion that the 
booms are passing through the wake of the rocket: in one revolution each of 

t 

the two booms passes through the wake once. It is assumed that this spin 
noise is an added component to the fine structure signal. Modelling the 
noise as band limited and consisting only of harmonics of the spin frequency 
opens the possibility of filtering out the spin noise, 

A second source of noise is the distortion added by the logarithmic 
electrometer. In the analysis leading to equations (4.1) and (4.2) it is 
assumed that £N e /<N e > « 1. Smith and Royrvik (1985) discusses the 
generation of harmonics by the electrometer for a sinusoidal input. He 
finds that the strength of each harmonic falls off rapidly as the order of 
the harmonic increases. The signal expected from the probe contains a large 
range of frequencies. If the power spectrum obeys a power law. k n , and 
has a large spectral index,, -n, then harmonics from the lower frequencies 
could be comparable to the strength of the spectrum at high frequencies. 

This would lead to overestimation of the spectrum at high frequencies. 

Smith and Royrvik (1985) concludes that for a power spectral index of -7 or 
greater the distortion is negligible. Spectra with smaller spectral indices 
are not expected from the fine structure experiment. However, spectra with 
exponential decreases as described in Section 3.5 could indeed be affected 
by this distortion. Fortunately, most of the information extracted from the 
spectra resides in the wave number ranges below the exponential decrease. 
Efforts to obtain information from the higher frequency ranges of a steeply 
sloping spectrum must not overlook the possibility of distortion. 

A third source of noise is interference from other experiments on board 
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the rocket. The fine structure experiment is especially susceptible to 
voltage transients that temporarily change the potential between the rocket 
body and probe. Other probe experiments on the rocket which use the rocket 
body as an electrode can generate such transients. Transients can also be 
generated by mechanical vibrations of the rocket. As an example, Table 4.2 
gives various events during a flight that could lead to electrical noise. 
Door ejection and boom erection refer to the insertion of the boom-mounted 
Langmuir probes into the ionosphere. These events are quite violent and 
vibrations can persist for about a second. An aural replay of the fine 
structure signal exhibits surprisingly distinct metallic vibrations 
corresponding to these events. Fortunately, the noise from these various 
sources of interference is easily identified, and although it is much 
stronger than the first two sources of noise described, its transient nature 
only causes isolated data dropouts as opposed to persistent contamination 
of the signal. 

4.4 Analog Prenrocessine 

An analog processing method developed by Klaus and Smith (1978) allows 
for easy identification of irregularities and a rough estimation of their 
spectra. The basic purpose of the system is to plot in various frequency 
bands the rms amplitude of the fine structure signal as a function of time. 

A block diagram is shown in Figure 4.4. The FM tape recorder plays back the 
combiner track, and the tunable discriminator decodes the desired channel. 
The 100 kHz reference signal is passed through a discriminator with a center 
frequency at 100 kHz. Tape speed fluctuations cause fluctuations in both 
the data signal and the recorded reference signal. The 100 kHz 
discriminator produces a voltage proportional to the tape speed fluctuation. 
Adding this difference signal to the discriminator output prevents 


Table 4.2 Flight log from 28 May Peru flight 


Nike Apache 14.532 


Site : 

Date : 
Launch Time 


Chi lea Peru 

28 May 1975 (Day 148) 

20 26 Universal Time (nominal) 
15 26 Local Standard Time 
20 25 s 59.7 (actual) 


Second Stage Ignition 
40 kft Baroawitch 
70 kft Barosvitch 
Door Ejection (61 km) 
Boom Erection (62 km) 
Probe Sweep On (128 km) 
70 kft Baroawitch Off 
40 kft Baroawitch Off 
Loas of Signal (Splaah) 


20 

26 

20 

26 

20 

26 

20 

26 

20 

26 

20 

27 

20 

32 

20 

32 

20 

33 


19.5 UT 

21.0 UT 

27.0 UT 

50.6 UT 

51.1 UT 

41.1 UT 

50.1 UT 
54.5 UT 

07.7 UT 



Figure 4.4 Analog data reduction system (Klaus and Smith, 1978) 
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fluctuations in tape speed from distorting the amplitude. The delay of the 
100 kHz discriminator is matched by a 75 psec delay added to the data line, 
preserving synchronicity between the tape speed compensation and data 
signals. Due to the FM nature of the tape, jitter in the tape speed causes 
jitter in both amplitude and phase of the output. The phase jitter is 
uncorrectable by analog means and must simply be tolerated as a source of 
noise. The compensated signal is then fed to a bandpass filter. The 
frequency bandwidths are chosen to increase geometrically (e.g., 10-20 Hz, 
20-40 Hz, 40-80 Hz, etc.). This gives spectral content at evenly spaced 
points on a logarithmic scale. The bandpass filtered signal is fed to an 
ac-to-dc converter and the resulting rms voltage is recorded on a chart 
recorder. A time code is simultaneously plotted. The result is a plot such 
as that in Figure 4,5, from the 28 May 1975 Peru flight. The frequency 
bands are converted to wave number bands by multiplying the frequencies by 
2ir/v where v is the rocket speed. 

This plot gives a broad view of the data, allowing for quick compari- 
sons between structures at different altitudes. The regularly spaced spikes 
seen clearly in the 0.18 - 0.36 rad/m band result from another probe 
measuring plasma resonance. The spike at 70.3 km is attributed to 
interference from an unknown source. Magnification of the raw fine structure 
data at 70.3 km is shown in Figure 4.6. A similar extraordinarily large 
signal swing occurs at 84.8 km. The apparent extended bandwidth of this 
spike observed in Figure 4.5 arises from clipping of the signal. Since the 
28 May 1975 Peru flight is an equatorial flight the beginning of the 
electrojet is visible at 90 km. The regions of most interest are the 
altitudes between 75 and 80 km and the irregularity at 87.5 km. Once the 
regions of interest have been identified from this broad overview, detailed 



Figure 4.5 Fine structure signal filtered into w-ve number bands 










inspection may begin at the proper altitudes. 

4.5 Digitizer 

To perform further analysis the analog data must be digitized. Prior 
to the development of the present digitizer system all digitization was 
carried out at the NASA facility at Wallops Island, Virginia. Since the 
turnaround time for the digitization varied from weeks to months, a 
digitizer system based on an Apple II microcomputer vsb developed at the 
University of Illinois capable of an 8-bit resolution, a digitizing rate of 
5 kHz, end s continuous record size of 25 kilobytes. This system solves a 
number of problems with the old systan. The Wallops Island system returns a 
digitial tape in a format compatible with a PDP computer. To use these data 
on the CDC mainframe this tape must be converted to CDC format with a 
conversion progran (Mclnerney and Smith, 1984). The program is quite lengthy 
and complicated, requiring the user to have a fair degree of familiarity with 
the CDC Cyber tape systan at the University of Illinois. To obtain data from 
the converted tapes, another Cyber progran must be run to control the 
mounting and playback of the tapes. The data generated are in a rather bulky 
format, making disk storage expensive. 

In the new systan the turnaround time for digitization is of the order 
of days, and the data are stored conveniently on floppy disks rather than on 
tapes. The data are transferred to the Cyber via a telephone modem. 
Essentially, the new system shifts control of the data to the user, avoiding 
complications and delays. Familiarity with the analog tape system is 

necessary to the extent needed for the analog processing described in 

> 

Section 4.3. It is not suggested that the new system is suitable for all 
applications; it is most efficient when handling small amounts of data. The 
time overhead needed for the old system cart pay off in studies requiring 
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large quantities of continuous data. In smaller projects this overhead may 
dominate the expended effort, making the microcomputer system a more 
attractive choice. 

A block diagram of the digitizer system is shown in Figure 4.7. The 
recorder plays back the FM combiner tr ck, the 100 kHz reference, and the 
two time codes. The EMR tunable discriminator selects the appropriate data 
channel and corrects for tape speed fluctuations as described in Section 
4.4. In addition, the discriminator has a built-in low pass filter 
necessary for reducing the cross talk from adjacent data channels. As 
suggested by Table 4.1, this filter is set to 1050 Hz for HLIG channel 18. 
The tape speed compensated data are sent to the A/D card for processing by 
the Apple computer. To assist in sampling data at regular intervals with 
respect to the tape, the 100 kHz reference signal is converted to a logic 
signal by a phase locked loop circuit. The Apple repeatedly counts a fixed 
number of these logic pulaes and takes a sample, resulting in a constant 
sampling rate. 

Use of the Apple internal clock for timing is avoided for a number of 
reasons. First, fluctuations in tape speed would cause the data to be 
sampled at irregular intervale with respect to the tape time frame. The 100 
kHz reference exhibits the same fluctuations as the tape, providing a clock 
synchronized to the data. Second, the Apple 1.023 MHz clock is not regular: 
every 65th ysec an approximately 140 nsec delay is added for videoscreen 
control. Cifntiauing with Figure 4.7, to ensure that the digitization starts 
at a precise time, a 1 pulse per second (1 pps) clock signal is obtained 
from the Datum time decoder, giving a positive transition exactly on the 
second. This transition is used to trigger the digitization at a precise 
and reproducible instant. 
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The A/D card is made by Mountain Computer and ie capable of taking a 
sample every 9 ysec, The card can sample data on any one of 16 input 
channels. The input lines are bundled in a standard 26 wire ribbon cable. 

To minimize high frequency noise on this unshielded cable, the cable is 
wrapped in copper tape. The analog to digital converter on the card accepts 
an input in the range -5 to 5 V, converting this to a digital value in the 
range 0 to 255. The -5 to 5 V range is not adjustable, so any scaling of 
the input must take place external to the card. Software controlling the 
card is straightforward. Each of the 16 inputs is treated as a memory 
location in the Apple peripheral card I/O space (addresses $C080 - $C0F.F) . 

If the card is in slot S, channel X may be read from address $C080 + S * $10 
+ X. Each time the address is queried the card begins another sample. Two 
read commands are thus required to take a valid sample: one to start 
sampling and one to load the result into the Apple. 

Interfacing of the timing signals to the Apple is performed by a 6522 
Versatile Interface Adapter (VIA) on a peripheral card made by John Bell 
Engineering. A system diagram of the VIA is shown in Figure A. 8. The card 
occupies 16 bytes of address space in the Apple II peripheral card ROM space 
($C100 - $C7FF), The various functions on the 6522 can request servicing 
from the Apple by setting a bit in the interrupt flag register (IFR). Each 
bit in the IFR has a corresponding bit in the interrupt enable register 
(IER). Preconfigured high, this corresponding bit enables a set bit in the 
IFR to pull the Apple IRQ line low, initiating an interrupt. 

The first timing signal required by the digitizer is a trigger to begin 
taking data. This is derived from the 1 pps signal from the time decoder, 
attached to line CA1 of the VIA. By configuring the peripheral control 
register (PCR) bit 0 high any positive transition on GAl sets IFR bit 1. 
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Figure 4.8 6522 system diagram. 








IER bit 1 ie low, disabling IFR bit 1 from pulling IRQ low. The Apple, 
waiting in a loop which continually examinee the IFR, exits the loop when 
IFR bit 1 goes high, starting the sampling. 

The second timing signal needed by the digitizer is an interrupt every 
20 cycles of the 100 kHz reference (for a 5kHz digitizing rate). There are 
two timers on the VIA: TIMERl and TIMER2. TIMER2 is used to count the 100 
kHz reference pulses from the phase locked loop circuit, the output of which 
is attached to line PB6 of the VIA. Each negative transition causes the 
counter value in TIMER2 to decrement by one. VThen the counter value 
decranents from $0000 to $FFFF , TIMER2 sets IFR bit 5. Since IER bit 5 ie 
set high, the VIA pulls IRQ low. The starting values of the counting 
registers are written by the interrupt service routine to the TIMER2 low 
latch and high counter byte. The second act, writing to the high counter 
byte, causes the VIA to place the low latch value into the low counter byte, 
to clear IFR bit 5, and to start TIMER2 counting pulses. Appropriately 
refilling the counter every interrupt cycle results in an interrupt request 
every 20 reference cycles. 

A flow chart of the digitizer software is shown in Figure 4.9. The 
systen consists of a BASIC program named DIGITIZER CONTROL and a collection 
of assembly language subroutines called DIGIAUTOSTART.OBJO (complete listing 
in Appendix II.l). DIGITIZER CONTROL, listed in Figure 4.10, uses POKE 
statements to pass various parameters to DIGIAUTOSTART.OBJO. These include 
the number of samples to be taken and various addresses. DIGIAUTOSTART.OBJO 
is assumed to begin at $8000. However, a NOP instruction is placed at 
$7FFF, and line 470 calls this address. This odd starting point is due to 
an idiosyncrasy of the Apple II that prevents a BASIC program from calling 
address $8000. The first subroutine called by DIGITIZER CONTROL (called in 
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50 REM DIGITIZER CONTROL 

60 REM THIS PROGRAM CONTROLS THE ASSEMBLER PROGRAM DIGIAUTOSTART.OBJO . 

100 REM TO SAVE MEMORY SPACE, REMOVE ALL REMARKS BEFORE LOADING. 

200 D$= ,M ’ : REM D$ EQUALS CTRL D. 

400 POKE -32754,47: REM CHANNEL PARAMETER - A/D SLOT * 16 + A/D INPUT 
CHANNEL. 

410 POKE -32753,16: REM HIGH BYTE OF BUFFER ADDRESS. LOW BYTE ASSUMED ZERO. 

420 POKE -32752,168: REM LOW BYTE OF NUMBER OF SAMPLES TO BE TAKEN. 

430 POKE -32751,97: REM HIGH BYTE OF NUMBER OF SAMPLES. 

460 POKE -32748,128: REM HIGH BYTE OF INTERRUPT SERVICE ROUTINE. 

465 POKE 32767,234: REM PUTS A NOP AT $7FFF. 

470 CALL 32767: REM CALIR DIGIAUTOSTART.OBJO AT $7FFF. FIRST 

472 REM ROUTINE INITIALIZES THE BUFFER AND PERIPHERAL CARDS. 

480 CALL -32757: REM CALLS THE WAIT LOOP FROM WHICH INTERRUPT OCCURS. 

500 PRINT D$;"BSAVa R##WTYP@XXS,A$1000,L25000" 

600 END 


Figure 4.10 Integer BASIC program DIGITIZER CONTROL. 
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line 470) configures the VIA and creates a buffer to hold the sampled data. 
Before returning to DIGITIZER CONTROL , the subroutine waits for the positive 
transition on line CA1. Once triggered the subroutine starts the counter on 
the VIA and returns to DIGITIZER CONTROL. The next subroutine called is 
simply a wait loop (called in line 480). When the VIA counter finishes 
counting down, IFR bit 5 is set, pulling IRQ low. An interrupt service 
routine takes a sample from the A/D card, restarts TIMERS, places the sample 
in the buffer, and returns to the wait loop. Since five seconds of data are 
needed, the buffer is 25000 ($61A8) bytes long. When the final sample is 
taken the interrupt service routine does not return to the wait loop but 
rather to DIGITIZER CONTROL. The buffer is then saved on floppy disk in 
binary format. 

Once the external hardware of the digitizer system is set up, running 
the program only requires loading the appropriate files into the Apple and 
starting the sampling at the desired time. The two files needed are 
DIGITIZER CONTROL and DIGIAUTOSTART.OBJO . Omission of address parameters in 
loading DIGIAUTOSTART.OBJO causes the file to be loaded at the default 
address $8000. Loading at another address necessitates changing the POKE 
statements of DIGITIZER CONTROL. An alternate starting address also 
requires the re^isembiy of the assembly code of DIGIAUTOSTART.OBJO. A 
manor y map showing the addresses used by the digitizer is given in Figure 
4.11. To take a 5 second data segment at, for example, 50 seconds into a 
flight, the analog tape is started about 20 seconds earlier. When the time 
code translator displays 49 seconds DIGITIZER CONTROL is run. Digitization 
then starts at exactly 50 seconds, the next positive transition on CA1 . 

After five seconds the sampling stops and the buffer is written to disk. 

The file name has the following format: "R" + rocket number + file type + 
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Figure 4.11 Memory map for digitizer software 



71 


M @’' ♦ elapsed flight time from launch to first sample + "S". For example, a 
data file for the 28 May 1975 Peru flight (code number 14.532) at 50 seconds 
into the flight would be labelled R14532DAT@50S. This name must be 
substituted into DIGITIZER CONTROL before the program is run. 

To achieve the highest resolution possible the data signal is amplified 
to the threshold of saturating the V input range. This amplification 
is provided by the discriminator. To calibrate the amplified data the 
digitizer is run with a frequency generator substituted for the FM data from 
the tape recorder. A block diagram is given in Figure 4.12. Since the 
frequency generator does not suffer from tape speed fluctuations, the 100 kHz 
reference is disconnected from the EMR discriminator and the delay is 
omitted, eliminating tape speed compensation. The frequency generator 
mimicks the output of a VCO. Stepping the frequency generator from 7.5% 
below to 7.5% above the center frequency of the VCO simulates stepping the 
VCO input from 0 to 5 V. The digitized values for a number of frequencies 
are taken with the A/D system, treating the generator as a data signal. 

This provides a calibration between VCO input and digitized output . Since 
the frequency generator does not exhibit tape speed fluctuations the tape 
Bpeed compensation is disabled. 

4.6 Data Postprocessing 

The final process in the data analysis is transferring the binary files 
from floppy disk to the Cyber mainframe computer and performing spectral 
analysis. Properly transferring and formatting the data iB a three-step 
process. The binary files must first be converted to text files. Next, 
these files are transferred to the Cyber via a 300 baud modem. Finally, on 
the Cyber the dafc* aust again be reformatted to simulate files created from 
Wallops Island tapes. Once the files are in the proper format they may ??e 
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Figure 4.12 Calibration for digitizer system 
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submitted to SPECTRA and other programs. SPECTRA estimates the power 
spectrum of an irregularity by processing the discrete Fourier transform 
(DFT) of a number of data segments. The estimate is then compensated for 
the filtering imposed by the experimental apparatus. The corrected estimate 
is plotted on a high resolution graphics terminal. The spectrum may be 
stored while others are calculated. Further processing can calculate 
averages, differences, and smoothing of the stored spectra. In an attempt 
to reduce the effects of spin noise, a simple moving average can be applied 
to the data before spectral analysis. In some cases this allows for a more 
confident estimation of the spectrum at lower wave numbers. 

To transfer a binary data file to the Cyber it must first be converted 
to a text file. The binary file is stored in memory starting at $1000. An 
APPLESOFT program called TEifTlFf (listing in Appendix II. 2), stored at $800, 
converts segments of 20 bytes infci» text strings 60 characters long. Single 
and double digit numbers are padded with zeros. These strings are collected 
into a text file and stored on floppy disk. Due to Apple DOS limits on file 
sizes, one binary file, 25000 bytes or 99 disk sectors, must be converted to 
two text files with a total of 302 sectors. 

These text files are then transferred to the Cyber via a telephone 
modem. The transfer is supervised by a set of programs called APPLE TERM 
(Version 3.0), written by Neil Homy for the Computing Services Office at the 
University of Illinois. The modern used is a 300 baud Micromodem II made by 
D.C. Hayes and Associates, Inc. It may be inserted into any of the 
peripheral card slots in the Apple. Booting APPLE TERM loads various 
subroutines and asks the user for the telephone number of the Cyber hookup. 
Once connected, the Apple acts as a Cyber terminal. The user opens a 
destination file for the text files using the Cyber's ICE editor. A useful 
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file name format ia "D" + the last three digits of the flight number ♦ "A" + 
the elapsed seconds from launch corresponding to the first datum in the file 
(e.g., D532A50). After opening the file the insert mode is invoked. The 
Cyber displays a line number and a question mark, the insert mode prompt. 

To start the transfer the user types <ESC> T. APPLE TERM then takes control 
of the keyboard, opens the specified text file, and sends the data. Each 
line of the ICE file receives one record of the text file. Upon completion 
of the transfer the user eaves the ICE file. 

At this point the ICE file is still not in the format acceptable to 
many programs developed for rocket data analysis. These accept the Wallops 
Island format. In this format, diagramed in Figure 4.13, the data are 
segmented into records, each with 1231 bytes. Each record consists of one 
real byte giving the starting time of the record, followed by 1250 integers 
divided into 250 frames. A frame has five channels, each of which is 
assigned to a different experiment. This channel arrangement facilitates 
synchronous comparisons between two or more experiments. Thus each record 
contains 250 samples from a single experiment. If data from only a single 
experiment are needed, this format is rather uneconomical in memory usage 
when compared to the ICE files. 

Conversion of the ICE files to the Wallops Island format requires a 
simple FORTRAN program called APL2CDC (complete listing in Appendix II. 3). 
Prograa APL2CDC reads in data from the ICE file, inserts the proper time 
bytes, packs each data point into a frame with four filler points, and 
writes out the constructed records to a new file. The unformatted write 
statement used writes out a whole 1251 point array at once. Other programs 
reading this file must use a corresponding unformatted read statement with 
an identically sized input array. Files created with unformatted write 
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Figure 4.13 Format of a data record using the Wallops Island format. 
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star/XDents cannot be edited by ICE. Due to the approximately fivefold 
increase in file size* limits on the Cyber system force the new file to be a 
direct access file, as opposed to the original indirect access file. The 
only major difference between the two file types li^s in the way the Cyber 
stores them. In a FORTRAN program, both file types are read from and 
written to in the same manner. Since the Wallops Island format entails an 
appreciable increase in storage costs and required memory, long-term data 
storage is most efficient using the ICE files. 

Once in the appropriate format th& data may be Fourier analyzed and the 
power spectrum estimated. This iB accomplished with an interactive FORTRAN 
pro gran called SPECTRA (complete listing in Appendix II. 4). A flow chart of 
S!*?ECTRA is shown in Figures 4.14 and 4.15. SPECTRA, written by Bruce Toraei 
(in Roth, 1982b), employs a power spectrum estimation method developed by 
Welch (1967). Ir. this method consecutive, possibly overlapping, 
periodograns are averaged together. An introduction to various methods of 
power spectrum estimation is available in Oppenheim and Schafer (19/5). 
SPECTRA requires access to two software libraries available on the Cyber: 
International Math and Statistical Library (IMSL) and the Graphics 
Compatibility System (GCS) library of plotting routines. SPECTRA first asks 
a number of questions allowing the user to control the analysis. These 
include the starting time of the data of interest, the number of data points 
desired, and the DFT block length. Tbp. starting time is approximately known 
from the analog preprocessing. A plot of the digital data allows for an 
even more precise determination of a starting point. Such plots are made by 
a FORTRAN program called DATPLOT (complete listing in Appei .x 11.5). The 
number of points to which SPECTRA is applied is also determined from a 
DATPLOT output. Figure 4.3 is an example of a plot by DATPLOT. The time 
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Figure £.15 Flow chart for program SPECTRA 
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span of interest is multiplied by the digitizing rate to give the number of 
points to be analyzed. This aequence of points may be divided into smaller, 
possibly overlapping segments. It is these segments to which the Fourier 
analysis is applied. 

The chosen length of the DFT may either match or exceed an individual 
data segment length, an inadequately long data segment being padded with 
zeros. When the data are segmented, the digital values are calibrated to 
give values of AN /<N >. A Hanning window is applied to each segment to 
produce a better power spectrum estimate. The DFT of each windowed data 
segment is found with a fast Fourier transform (FFT) routine in XMSL. The 
squared magnitude of ths DFT times some scaling constants gives the power 
spectrum estimate. To correct for the various filters in the rocket 
experiment, SFECTBA divides the spectrum by the theoretical transfer 
function of the filters. This corrected spectrum is plotted using GCS 
graphics routines. The individual spectra are then averaged together to 
give a best estimate of the power spectrum. The averaged spectrum may be 
stored while a different averaged spectrum is calculated. Two spectra may 
be stored at one time. 

A number of modifications to SPECTRA allow for further processing of 
the spectra. The aim of these processes is to produce a better estimate of 
the spectral index. First, two stored spectra may be averaged together. 
This allows for the spectra of two nunconsecutive segments to be averaged. 
Such an average may be necessary to avoid a data segment flawed by inter- 
ference. Second, the difference of two spectra may be calculated. With 
this function an estimate of the background noise power spectrum may be 
subtracted from another spectrum with both signal and noise. Third, a 
spectrum may be smoothed by dividing the frequencies into groups whose size 
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increases with frequency by a constant factor. This has the effect of 
dividing a logarithmic frequency axiB into equal segments. The correspond- 
ing magnitudes in each frequency range are averaged and plotted. 

To study the effects of spin on the power spectrum, a moving average is 
applied to the raw data prior to spectral analysis. This process is per- 
formed by a FORTRAN program called SPINFIX (complete listing in Appendix 
11.6) which convolves an entire file of data points with a rectangular 
window function. The window has height of unity and length equal to the 
spin period. In the frequency domain this process transforms to multiplying 
the power spectrum by T^ sin (o>T/2 ) /(u>T/2 ) , where T^ is the, number 
of data points in a spin period, T is the spin period in seconds, and u is 
the angular frequency. The filtering action of the averaging necessitates 

some minor changes to SPECTRA. Multiplying the calculated power spectrum by 

2 2 
( Tif T/T ) leaves a spectrum multiplied only by sin (irfT) where f is the 
P 

frequency. ThiB should eliminate strong peaks at multiples of 1/T, the spin 
frequency. The power spectrum of the data in Figure 4.3 is shown in Figure 
4.16. Figure 4.17 shows the power spectrum of the same data after SPINFIX 
is applied. The FFT length is chosen to extend over at least one spin 
period. Figure 4.17 shows a dramatic reduction of spectral power in the 
lower frequency ranges. 

4.7 Summary 

The analysis described in this chapter may be viewed as a series of 
refinements toward the ultimate goal of estimating the power spectrum and 
the spectral index. Irregularities must first be identified and distin- 
guished from noise added by the experimental apparatus. The identification 
is achieved by the analog processing system developed by Klaus and Smith 
(1978). This system also provides a rough estimate of the fine structure 
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spectrum. Once interesting portions of the signal are identified, the data 
may be digitized. Detailed plots of the digital data provide the exact 
locations of irregularities within the data record. These irregularities 
are then subjected to digital analysis. The resulting spectrum has far 
better resolution than that afforded by the analog analysis. Smoothing and 
noise subtraction provide further ways of improving the spectral index 
estimates. 
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5. RESULTS 


5.1 Introduction 

This chapter presents some results of the analysis described in Chapter 
4. Fine structure data are examined from the nine daytime rocket flights 
since 1975 that have carried the fine structure experiment. Possible 
ionospheric irregularities are found with five of the nine flights, and 
power spectra of the irregularities are calculated. Many spectra exhibit a 
wave number range with a spectral index close to -5/3, consistent with the 
existence of an inertial range of turbulence. The slope and wave number 
range of the best fit lines to the spectra are tabulated. Of particular 
interest is the upper wave number extent of the best fit line. This limit 
is compared to predictions made by Raetogi and Bowhill (1976). 

Section 5.2 describes each rocket flight on which irregularities are 
observed. A particular irregularity from the 28 May 1975 Peru flight iB 
analyzed to illustrate the analysi •< procedure. Measured parameters of all 
observed irregularities are tabulated. Section 5.3 discusses the results 
derived in Section 5.2 and makes comparisons to the spectral model of Hill 
and Bowhill (1976) and to the inertial range predictions of Raetogi and 
Bowhill (1976). Section 5.4 presents the conclusion that the observed 
irregularities are consistent with those expected from atmospheric 
turbulence. 

5.2 Observations of Irregularities 

The 28 May 1975 Peru flight was launched from the Punta Lobos rocket 
range at Chilca., Peru at 15 26 local standard time (LST). The coordinates of 
the site are 12.5° S and 76.8° W. The prime objective of the flight «as 
to observe the equatorial electrojet. This site is particularly interesting 
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since the Jicamarca radio observatory is located only 60 km north of the 
rocket range. The close proximity of the two sites facilitates coincident 
rocket and radar measurements. More details of the launch operations and 
flight are available in Klaus and Smith (1978), including comparison of 
Langmuir probe electron density measurements to incoherent scatter 
measurements taken at Jicamarca. 

The fine structure data from the 28 May 1975 Peru flight, although 
containing some interference, in general are of exceptional quality. There 
is no discernible spin contamination since a nose tip probe was used. 

Analog preprocessing produces the strip chart records shown in Figure 4.5. 
The data are marred by interference spikes appearing every V»!f second, 
attributable to an RF resonance probe on the rocket. Two isolated noise 
spikes occur at 70.3 and 84.8 km. Their source is unknown, but the 
shortness and strength of the spikes suggest they are due to interference. 
Not so easily explained are the four fluctuations seen following the large 
spike at 84.8 km. Examination of the digital record reveals no obvious 
common features between them. They cannot be attributed to Bpin since their 
spacing is not exactly regular or close to the spin period. Aside from 
these three groups of isolated noise spikes, the only other noise present is 
the unavoidable noise due to the telemetry system. 

The fine structure data reveal the presence of six irregularities 
between 70 and 90 km. Table 5.1 lists their altitudes and vertical extent. 
To illustrate the analysis process, the irregularity at 74.9 km is 
considered here in detail. The top plot in Figure 5.1 shows a portion of 
the digital record at 74.9 km. This represents the raw fine structure 
signal without correction for the various filters in the experiment. In 
particular the observed signal is not corrected for the low frequency 


Table 5.1 Parameters of irregularities for nine rocket shots. No 
irregularities are found with the last four shots listed 
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Figure 5.1 Upper plot shows portion of digital record from the 

28 May 1975 Peru flight. Launch time at 15 26 (LST). 
Altitude of the irregularity center is 74.9 km. Lower 
plot shows the magnetometer output. 
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attenuation from the ac amplifier. Low frequency components thus appear 
de-eaphasized. The upper independent axis gives relative distance along the 
flight path while the lower axis shows elapsed time from launch. The center 
of the irregularity, at 900 m on the distance scale, occurs at 74.9 km 
altitude. Interpretation of the observed temporal fluctuations as spatial 
fluctuations with scale sizes read directly from the distance axis requires 
the assumption of isotropy. With this assumption, a temporal oscillation 
with frequency f corresponds to a spatial oscillation with wave number k 
given by 

k - 27tf/v (5.1) 

where v is the rocket velocity. Host spectra presented in this chapter, 
although calculated for temporal frequencies, display a wave number axis 
with values found from equation (5.1). 

The two-sided character of the !•+.? km irregularity suggested by Figure 
4.5 shows up quite strongly in Figure 5.1. It might be suspected that the 
two distinct portions of the irregularity are the result of rocket spin. 
However, analysis of the irregularity reveals two countering arguments. 
First, SB observed in Figure 4.3, each cycle of the spin modulation should 
closely resemble the neighboring cycles. This is not the case in the top 
portion of Figure 5.1. Second, any matching features of the two suspected 
cycles should be separated exactly by the spin period. A portion of the 
magnetometer signal coincident with the irregularity is shown in the bottom 
portion of Figure 5.1. The spin period is measured to be 0.245 seconds. In 
contrast, the time displacement between the two large dips on the left of 
each apparent spin cycle in Figure 5.1 is 0.207 seconds. 

Two power spectra are estimated, one for each of the two prominent 
halves of the irregularity. The two power spectra each have an FFT length 
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of 512 points. The average of the two spectra is shown in Figure 5.2a. The 
frequency scale is converted to a wave number scale using equation 5.1 with 
n speed of 1.50 x 10 J m/sec. The units on the vertical axis are converted 
from Hz*** to (rad/m)"*. A noise floor due to the telemetry system is 
seen at approximately 10** 7 (rad/m)~*. This limits the useful wave 
number range to below 1 rad/m, where the irregularity spectrum crosses the 
noise floor. To help measure the spectral index, a noise estimate, Figure 
5.2b, is subtracted from the averaged spectrum. The difference is shown in 
Figure 5.2c. The data used in calculating the noise estimate come from a 
data segment adjacent to the irregularity. The observed noise is attributed 
to the telemetry system. Since the telemetry noise enters the signal after 
the ac amplifier, the SPECTRA correction for the high pass filtering boosts 
the low frequency telemetry noise components above their actual level. This 
partially accounts for the spectral strength at low frequencies observed in 
Figure 5.2b. The difference spectrum is smoothed by SPECTRA ueing the 
process discussed in Chapter 4. The window multiplication factor is 1.4. 

The resulting smoothed spectrum is shown in Figure 5. 2d. 

Three separate regions of the smoothed spectrum may be distinguished. 

A line may be fitted between the first point on the left and the point at 
0.23 rad/m. The best fit line has a slope of -1.7. The second region, 
between 0.23 and 2.6 rad/m, exhibits a slope of *-3.7. The third region is 
that above 2.6 rad/m. In this region telanetry noise overwhelms the signal 
and little information may be extracted. The V-shaped dips shown in Figure 
5. 2d are an artifact of the computer plotting routine, which rounds up 
off-scale points to the lowest point on the graph. The average of a set of 
points on the difference spectrum (Figure 5.2c) may be negative and hence 
off-scale. Negative points are therefore set to the lowest point on the 
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Figure 5.2b Estimated noise power spectrum for the 74„9 km irregularity 
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Figure 5.2c Difference between the irregularity power spectrum and the noise 
power spectrum shown in Figures 5.3a and 5.3b, respectively. 
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vertical axis. 

Figure 5.3 showB the smoothed power Bpectum of an irregularity at 87.5 
km. The best fit line to the wave number region below 0.19 rad/m has a 
elope of -1.8. The second steeper slope has a value of -7.5. The total 
power is much greater than that in the irregularity at 74.7 kxu. Figure 5.4 
shows the electron density profile from the 28 May 1975 Peru flight. Note 
that the altitude of the 87.5 km irregularity coincides with a large 
electron density gradient. It might be suspected that the irregularity is 
merely the response of the ac amplifier to the sudden increase in electron 
density. However, the data record at times immediately following the 
irregularity shows very little activity, in spite of the continued presence 
of the gradient. It is thus reasonable to assume that the source of the 
irregularity is geophysical and not instrumental. 

The results described above are given in Table 5.1, along with the 
results from other irregularities from the 28 May 1975 Peru flight. The 
wave number range over which the power law extends is determined from the 
unsmoothed difference spectra. In most cases the best fit line extends to 
the lowest wave number on the spectrum. The low limit listed in Table 5.1 
thus reflects the lowest wave number for which the ijFT is calculated and is 
a function of the FF? size and the digitizing rate. The high limit listed 
is the largest wave number at which the spectrum remains on the best fit 
line. 

As mentioned previously, the digitized fine structure signal (e.g., 
Figure 5.1) is not corrected for the low frequency attenuation from the ac 
amplifier. To measure the maximum deviation of an irregularity, this 
attenuation must be accounted for. To correct the observed maximum 
deviations, the main frequency component of the deviation is estimated and 
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Figure 5.4 Electron concentration profile for the 28 May 1975 Peru flight 
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the attenuation calculated. For example, the maximum deviation observed in 
Figure 5.1 is 1.2%. The main frequency component of the deviation appears 
to be 44 Hz. With a lover 3 dB frequency of 54 Hz, the amplifier attenuates 
a 44 Hz signal by 1.6. The corrected deviation is thus 1.6 x 1.2% ■ 1.9%. 

All deviations in Table 5.1 represent corrected values. 

The 18 January 1976 Wallops Island flight was launched from Wallops 
Island, Virginia (37.8° N 75.5° W) at 14 12 LST. The primary purpose of 
the flight vae to observe the winter anomaly. The rocket also carried a 
mass spectrometer prepared by the University of Bern, Switzerland. The 
placement of the spectrometer beneath the nose cone prevented the use of a 
nose tip probe configuration. Instead, the probes were placed on booms 
extended from the sides of the rocket payload. Unfortunately, the nose cone 
shell detachment to expose the spectrometer malfunctioned during the flight, 
damaging the telemetry system and the propagation experiment. Thus no 
calibration between electron density and probe current could be calculated. 
The probe and fine structure experiments were cot damaged. A description 
of this flight is given in Edwards (1976) . 

Since no calibration is available for the probe current the electron 

density can be only roughly estimated. Processed values for the probe 

current are also unavailable. Fortunately, the raw probe current signal is 

available on analog tape (IRIG channel 15). Plotting this signal on a strip 

chart provides a way to measure the probe current, albeit with less accuracy 

than a digital analysis would yield. To convert the observed current to 

electron density, an average calibration value from other flights is applied. 

9 "3 

Using a calibration value of 3 x 10 cm /A yields the concentration 
data shown in Table 5.1. Due to the coarseness of the probe current 
measurement and the high degree of spin contamination, calculations of the 


98 


electron density grad .ent are unreliable, and only order of magnitude 
estimates are given. It is possible to observe a large positive gradient in 
the probe current beginning at 84 km. 

The data reveal four ionospheric irregularities. Their altitudes and 
thicknesses are given in Table 5.1. As observed in Figure 4.3, spin 
contamination is a serious problem in the fine structure data from the 18 
January 1976 Wallops Island flight. To calculate a power spectrum estimate, 
a number of noise samples are analyzed and their average spectrum subtracted 
from the irregularity spectrum. The segments of noise and irregularity data 
have the same phase with respect to the spin contamination. The slope and 
limits of the best lines fit to the spectra are given in Table 5.1. A 
smoothed spectrum from the irregularity at 82.2 km is given in Figure 5.5. 

The 12 August 1976 Wallops Island flight was launched at 10 54 (l£T). 
The main purpose of the probe experiment on the rocket was to observe 
sporadic E layers. Details of the probe current analysis are given in 
Edwards (1977) . The fine structure e::periment again used a boom probe 
configuration. The spin contamination was very small. It should be noted 
that the propagation experiment data on this flight are unreliable above 84 
km (Edwards, 1977). This unreliability contributes uncertainty to the 
electron density calculations above 84 km since the propagation experiment 
is used to calibrate the probe current. 

Only one irregularity is observed in the data, occurring at 83.5 km. A 
smoothed spectrum from this irregularity is shown in Figure 5.6. This 
irregularity is superposed upon the steep electron density gradient often 
observed between 80 and 90 km. The peak fluctuation of the irregularity is 
unusually high: 13%. 

The 8 January 1977 Wallops Island flight was launched at 14 03 (LST). 



Smoothed difference 
18 Jiinuary 1976 Wal 
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The primary purpose of the flight was to observe the winter anomaly. A 
description of the flight preparations is given in Edwards (1977) . The data 
suffered from moderate spin contamination. 

Only one irregularity was observed. As shown in Table 5 , 1 £ Hie 
irregularity extended over a relatively short altitude range. Figure 5.7 
gives the smoothed spectrum. The large dip at 0.06 rad/m results from the 
estimated noise spectrum exceeding the estimated signal plus noise spectrum. 

The 27 February 1983 Peru flight was launched from the Punta Lobos site 
at 11 33 (LST). Its principal purpose was to study electron concentration in 
the equatorial D and E regions. On this flight' the fine structure 
experiment was modified to allow greater gain in the ec amplifier. Those 
changes are documented in Roth (1982a, 1982b). The flight did not carry the 
mass spectrometer, permitting use of the nose tip probe configuration. 

The data reveal one irregularity at 85.8 km, superposed on a steep 
electron density gradient. This irregularity is thoroughly investigated by 
Royrvik and Smith (1984) and is compared to simultcneous coherent scatter 
measurements at Jicamarca. A spectrum of this irregularity is shown in 
Figure 5.8. 

5 .3 Discussion of Results 

All spectra calculated show a low wave number rang?, over which a power 
law may be fitted. That the slope of this fitted line usually has a value 
close to -5/3 (Table 5.1) suggests the presence of an inertial range of 
turbulence. As a consequence of spectral aliasing, the one-dimensional 
spectra do not show the spectral decrease expected below the inertial range 
in the three-dimensional spectra. In the wave number ranges above the 
suspected inertial range, the measured spectra do not always demonstrate the 
same behavior. Some show a second power law range, as in Figures 5.2a and 
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Figure 5.7 Smoothed difference spectrum of the 75.9 km irregularity from the 8 
January 1977 Wallops /eland flight. 
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5.3, while others show a more curved decrease in spectral power. In all 
spectra the telemetry noise floor prevents observation of the upper wave 
number behavior over any sizable wave number range. It is interesting to 
note thet three of the five flights discussed show irregularities superposed 
upon the electron density gradients. All three irregularities have quite 
high maximum AN e /<N e > values. 

Qualitative comparisons may be made with the spectral model of Hill and 
Bowhill (1976) described in Chapter 3. Figure 3.3 shows the predicted 
spectra for various values of the parameter 5, equal to 2£(k*/k c ) f/ . 

This parameter should be of the order of unity in the D region. The 
predicted spectrum for € - 1 shows no viscous-convective Inge or power law 
above t } .j inertial range. No viscous-convective range is found in any of 
the irregularities observed, in agreement with the model. Figures 5.5, 5.6, 
5.7, and 5.8 show a curved decrease as predicted by the spectral model. In 
contrast to the model. Figures 5.2a and 5.3 show a power lav dependence 
above the apparent inertial range. 

The results of Table 5.1 may be compared with predictions by Rastogi 

and Bowhill (1976) of the wave number extent of the -5/3 power law 

dependence. For two models of e, Rastogi and Bowhilx (1976) calculate the 

variation of k^ with altitude. Figures 5.9a and 5.9b, using model 1 and 

model 2 of e , respectively, show three wave number ranges. Below ky the 

—5/3 

energy spectrum should vary as k . Above k^ the spectrum is pre- 
dicted to fall off exponentially. Superposed upon Figures 5.9a and 5.9b are 
the upper wave number limits .given in the last column of Table 5.1. It 
should be noted that the calculated spectra used in creating Table 5.1 are 
one-dimensional spectra. As suggested by Figure 3.1, the upper limit of 
the -5/3 power law reaches a higher wave number in the three-dimensional 
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Figure 5.9a Suggested profile of cutoff wave number demarcating range 
where —5/3 power law holds (k < k v ). Wave number of 
Jicamarca (50 MHz) and Arecibo (430 MHz) incoherent 
scatter radars are shown. The wave number range used for 
VHF forward scattering is also shown, e from model 1 
(Rastogi and Bowhill, 1976). 
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Figure 5.9b Suggested profile of cutoff wave number demarcating range 
where “5/3 power law holds (k < k v ). Wave number of 
Jicamarca (50 MHz) and Arecibo (430 MHz) incoherent 
scatter radars are shown. The wave number range used for 
VHF forward scattering is also shown. e from model 2 
(Raetogi and Bowhill, 1976). 
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scalar spectrum than in tha one-dimensional spectrum. The values in Table 
5.1 thus underestimate the upper wave number limit o£ the apparent inertial 
range. Also, the values plotted in Figures 5.9a and 5.9b are from scalar 
spectra, as opposed to energy spectra. However, viscous and diffusive 
effects should be expected to begin at approximately the same wave number 
since k c - k d in the mesosphere. It is thus not surprising that the 
scalar spectrum inertial range extends to the same order of wave number as 
the energy spectrum as suggested by Figures 5.9a and 5.9b. 

5 .A Cone lusions 

The scalar spectra shown in this chapter demonstrate consistency with 

theories of homogeneous, isotropic turbulence measured in one dimension. 

All the spectra show a definite wave number range over which the scalar 

spectrum varies as a negative power of k. Most theories of turbulence 

- 5/3 

predict an inertial range with a k dependence. The data in Table 5.1 
give spectral slopes close to -5/3, consistent with an inertial range. 

Above the inertial range, turbulence theory predicts that viscous and 
diffusive effects cause a precipitous decrease in the spectrum. Due to the 
small Schmidt number, no viscous-convective range is expected in mesospheric 
turbulence. 

Although not in agreement on the exact form, most turbulence theories 
predict an exponential decrease in the spectrum above the inertial range. 
Hill and Bowhill (1976) incorporate these theories into a model for the 
scalar spectrum. The model is extended to one-dimensional spectra under the 
assumption of isotropy. Many of the measured spectra agree qualitatively 
with the shape of the predicted spectrum. To show further consistency with 
turbulence theory the spectra are compared to predictions by Rastogi and 
Bowhill (1976). They predict the upper wave number extent of the inertial 


range. The measured limitB agree approximately to the prediction*. Perhaps 
more importantly, the measured inertial limits decrease with increasing 
altitude, in agreement with the functional dependence on altitude predicted 
by Rastogi and Bowhill (1976). The agreement of the measured spectra with 
the models of Hill and Bowhill (1976) and Rastogi and Bowhill (1976) support 
the conclusion that the observed ionospheric irregularties are consistent 
with irregularities expected from atmospheric turbulence. 
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6. SUMMARY AND SUGGESTIONS FOR FUTURE WORK 


6 .1 Summary 

Mesospheric electron density irregularities are often observed by both 
radar and rocket studies. These irregularities are suspected to be the 
result of turbulent advection of ions and electrons. Possible mechanisms 
for generating atmospheric turbulence include Relvin-Helmholtz and Rayleigh-* 
Taylor instabilities. Turbulence may be qualitatively described as the 
cascading of energy from large scales to small scales. At the smallest 
scales the turbulent energy is dissipated by viscosity (for turbulence in a 
velocity field), or by molecular diffusivity (for turbulence in a scalar 
field). Most theoretical studies of turbulence attempt to estimate the 
spectrum of homogeneous, isotropic turbulence. Physically, such a spectrum 
gives the turbulent energy per wave number. Most theories predict an 
inertial range in the spectrum for flows with sufficiently high Reynolds 

_5/3 

number. The spectrum is expected to vary as k in the inertial range. 

The theories do not agree closely in describing the scalar spectrum above 
the inertial range, although most theories in general predict an exponential 
decrease with increasing wave number. 

Langmuir probes on sounding rockets provide aa in situ method of meas- 
uring mesospheric electron density irregularities. In the fine structure 
experiment flown by the University of Illinois, a Langmuir probe is extended 
from the rocket into the ionosphere. As the rocket passes through an irreg- 
ularity, the probe current fluctuates in time. With the assumption of iso- 
tropy, temporal fluctuations are interpreted as spatial fluctuations meas- 
ured in one dimension. The measured current from the probe is processed 
on board the rocket to produce a voltage proportional to AN /<N >. The 
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signal is telanetered to the ground and is stored on an FM analog tape. 

Using digital analysis, the power spectra of the irregularities are 
estimated. The analog data are digitized with a microcomputer based 
digitizer. The microcomputer stores the digital data on floppy disks. The 
data are then transferred to a CDC Cyber mainframe computer for spectral 
analysis. The FORTRAN program SPECTRA estimates the power spectra of the 
irregularities using an FFT routine available on the Cyber. Further 
processing subtracts estimates of the background noise spectra and smooths 
the result of the subtraction. These smoothed spectra help to identify 
spectral indices and general trends in the spectra. 

The measured spectra appear consistent with the predictions of turbu- 
lence theory. The observed spectra all exhibit power laws with spectral 
indices close to -5/3, consistent with an inertial range of turtulence. 
Comparisons to a spectral model by Hill and Bowhill (1976) and to predic- 
tions of inertial range limits by Rastogi and Bowhill (1976) show further 
consistency with theory. 

6.2 Suggestions for Future Work 

There are two significant ways by which the fine structure experiment 
could be improved. Coincident radar measurenents would help to establish 
the existence of a turbulent layer and perhaps help to determine isotropy. 
Reduction of telemetry noise would increase the upper wave number limit of 
the turbulence spectrum. 

Coincident radar measurements of an irregularity would provide data 
unattainable from the Langmuir probe experiment. The coherent-scatter radar 
technique can measure a single component of the scalar spectrum. The 
advantage of the radar technique is its ability to measure a spectral 
component that in the rocket data is obscured by telemetry noise. The radar 
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thus provides insight into the shape of the spectrum at wave numbers above 
the inertial range. The combination of radar and rocket measurements is 
used by Royrvik and Smith (1984) to compare measuronents of an irregularity 
from the 27 February 1983 Peru flight with coincident radar measurements at 
Jicamarca. Badar can also measure aspect sensitivity of au irregularity, an 
indication of anisotropy. This property cannot be measured by the fine 
structure experiment. Knowledge of isotropy is essential in allowing 
comparisons of measured one-dimensional spectra with theories predicting 
three-dimensional spectra. 

Reduction of telemetry noise would increase the highest observable wave 
number of a measured spectrum. Telemetry noise could be greatly reduced by 
use of a digital telemetry system. One possible analog nois& reduction 
scheme is described in Roth (1982b). Essentially, the method takes advantage 
of the relatively large spectral magnitude expected at low frequencies in 
an irregularity spectrum. The fine structure experiment described by Klaus 
and Smith (1978) uses an approximately 50 Hz 3 dB point in the ac amr/iif ier 
high pass filter to prevent saturation of the telemetry system by the low 
frequencies. Further attenuation of the low frequencies allows the overall 
gain to be increased, giving a better ratio of signal to telanetry noise. 

The new system places the 3 dB point at 400 Hz and raises the ac amplifier 
gain to 400. The power spectra calculated from this new systeu show an 
order of magnitude decrease in the telemetry noise floor. Ideally the 3 dB 
point should occur just below the frequency at which the irregularity 
spectrum crosses the telenetry noise floor. Inclusion of this system in 
future rocket flights would raise the highest observable wave number of the 
calculated spectra. 

The above analog noise reduction scheme could be improved by the use of 


a more elaborate filter. If instead of a single-pole low pass filter the ac 
amplifier used a filter more closely matched to the expected irregularity 
spectrum, then the amplifier gain could be increased further, improving the 
ratio of signal to telemetry noise. SPECTRA could then easily correct the 
calculated irregularity spectra for the filter transfer function. 
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APPENDIX I FREQUENCY RESPONSE OF AN AC AMPLIFIER 

To make corrections for the filtering of the ac amplifier used in the 
fine structure experiment the upper and lover 3 dB points of the frequency 
response need to be known. Figure 1.1 shows a circuit diagram of the 
amplifier circuit. The logarithmic electromecer driving the ac amplifier 
has an output impedence of 1 k% Including the electrometer output imped- 
ance in series with the 11 kft resistor gives an amplifier gain of 100 and 
upper and lower 3 dB frequencies of 1.7 kHz and 60 Hz, respectively. 

The frequency response of three amplifier circuits was measured. To 
simulate payload conditions the amplifiers were tested after their incor- 
poration into circuit modules used in an actual payload. A Hewlett Packard 
frequency synthesizer (model 3325A) was used as the input to the ac ampli- 
fiers. The frequencies at which the amplifier output was attenuated by 
- 1/2 

2 times the maximum output are given in Table 1.1. 

The deviation of the measured values from theoretical predictions does 
not seriously affect the data analysis. The measured lower 3 dB frequencies 
are in good agreement with the value of 54 Hz found by Klaus and Smith 
(1978) for the 28 May 1975 Peru flight. Unfortuuately , the lower 3 dB 
frequencies were not measured on the other eight flights studied in this 
thesis. The lower 3 dB frequency is needed by SPECTRA to correct for the ac 
amplifier high pass filter. The value 54 Hz is adopted to represent the 

3 dB point for all flights. Although the measured upper 3 dB frequencies 
differ markedly from the expected values, telemetry noise dominates the 
irregularity spectrum at frequencies of the order of the upper 3 dB fre- 
quency. Exact correction for the ac amplifier low pass filtering is thus 
unimportant. However, should telemetry noise be reduced on future flights 
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to such a level that the Irregularity spectrum is stronger than the high 
frequency telemetry noise, then pains should be taken to accurately measure 
the upper 3 dB frequency. 


APPENDIX II SOFTWARE LISTINGS 


II. 1 DIGIAUTOSTART 


SOURCE FILE: 
0000 1 
0000: 

0000 : 

OOQOt 

0000: 

0000: 

0000: 

0000 : 

0000: 

0000 : 

0000: 

0000 ; 

0000: 

0000: 

0000: 

0000: 

0000: 

0000: 

0000; 

0000; 

0000: 

0000: 


DIGIAUTOSTART 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 


THIS PROGRAM TAKES A LARGE NUMBER 

OP SAMPLES PROM THE A/D CARD AND PUTS THEM IN A BUFFER, THE STARTING 
ADDRESS is G 1 VQ) BY BUFH AND BUPL. SINCE BASIC MEMORY STARTS AT $9600 
(WITH DOB LOADED) AND WORKS DOWN, DIGIAUTOSTART SHOULD BE LOADED A WAYS 
BELOW THAT TO ALLOW POR A BASIC PROGRAM. $8000 IS SUGGESTED. THE PROGRAM 
FROCRAM DIGITIZER CONTROL TAKES UP ABOUT A PAGE OP MEMORY ($100). AS IS, 
THIS PROGRAM RUNS AT A DIGITIZING RATE OP 5 KHZ. THE JOHN BELL INTERFACE 
CARD IS IN SLOT 4. THE A/D BLOT IB NOT FIXED BUT IS CONTROLLED BY CHNL. 

THIS PROGRAM ALSO ASSUMES A REFERENCE FREQUENCY OR 100 KHZ, OR 10 USEC 
PER REFERENCE PULSE. TO CHANGE THE DIGITIZENG RATE THE OGPAUS ASSIGNMENT 
MUST BE CHANGED. TO CHANGE THE J.B. CARD BLOT THE BELBLT ASSIGNMENT MUST 
BE CHANCED. BOTH CHANGES REQUIRE RECOMPILING. THE DATA COLLECTION IS 
IS TRIGGERED BY A POSITIVE TRANSITION ON THE CAl LINE OF THE JOHN BELL CARD. 
A 1 PUISB PER SECOND SIGNAL IS AVAILABLE FROM THE TIME CODE TRANSLATOR 
IN THE ROCKET LAB. THE PROGRAM USES THE JOHN BELL CARD TO COUNT A. 

NUMBER OF EXTERNAL CLOCK CYCLES. WHEN THE COUNT IS UP, THE J.B. CARD 
GENERATES AN AN INTERRUPT REQUE3T SIGNAL. THE COMPUTER GETS OUT OF THE 
WAIT LOOP AND JUHPS TO THE INTERRUPT SERVICE ROUTINE CALLED IHNDLR. 

THIS PROGRAM IS A MODIFICATION OP A PROGRAM CALLED CLKREAD. CLKREAD 
AND OTHER A/D PROGRAMS CAME WITH THE MH'WIAIN HARDWARE D/A A/D CARD, 
DIGIAUTOSTART BY ROD ST0LTZPU8 , WRITTOJ SEPTEMBER 1984 
FOR THE AERONOMY LAB AT THE UNIVERSITY OP ILLINOIS. 


NEXT OBJECT FILE NAME IS DIGIAUTOSTART. OB JO 


8000: 


23 

ORC 

$8000 

jPUT PROGRAM AT LOCATION $8000 

8000:78 


24 START 

SEI 


■.CALLED FROM A BASIC PROGRAM, START JUMPS TO 

8001: 


25 ; 



SUBROUTINE IH1T AND INITIALIZES THE PERIPHERAL CARDS 

8001 :4C 72 90 

26 

JMP 

IN IT 


8004:76 


27 STOP 

SEI 


-.RESTORES REGISTERS AND SUCH. 

8005 :4C El 

80 

26 

JMP 

FINISH 


8008 s4C IB 80 

29 TICK 

JMP 

IHNDJR 

jINTERHUPT SERVICE ROUTINE. THE INTERRUPT 

800B: 


30 ; 



VECTOR $03FE AND $03?P CONTAINS IHNDLR' S 

BOOB: 


31 ; 



ADDRESS NOT TICK'S ADDRESS. 

800B-.4C 07 

81 

32 DUMMY 

JMP 

WAIT 

jA LOOP IN WHICH TO WAIT FOR AN INTERRUPT. 

OOOE: 


33 ;**»****>*« **************** ***** *** * ***** «»******>****<! ** ******* 

8003: 


34 ;DATA STORAGE LOCATIONS: POKED BY A BASIC PROGRAM. 

800E; 


35 .*****<m***«************************«*******«**««*************** 

800E: 


36 CHNL 

DS 

1 

{■A/D SLOT# * 16 + A/D CHANNEL NUMBER (0 TO 15). 

3001? : 


37 BUFFER 

DS 

1 

{HIGH BYTE OF STARTING ADDRESS OF THE BUFFER 

8 D 10 : 


38 ; 



WHICH HOLDS THE DATA. LOW BYTE ASSUMED TO BE $00. 

8010: 


39 ; 



THIS ASSUMPTION MAY T-, ALTERED BY CHANGING BUFL 

8010: 


40 ; 



IN THE INIT SUBROUTINE. 

8010: 


41 COUNT 

DS 

2 

{NUMBER OF SAMPLES TO BE TAKEN AND STORED. 

8012:EA 


42 

NOP 



8013:EA 


43 

NOP 



8014: 


44 INTADH 

DS 

1 

{HIGH BYTE OF IHNDLR ADDR. PROBABLY EQUAL TO THE 

8015; 


45 ! 



HIGH BYTE OF THE ORG ADDRES 8 (SEE FIRST LINE). 

8015: 


45 .A************* ************************************************* 

8015: 


47 ; INTERNAL VARIABLES 

NEXT 

8015: 


(|0 JWWWWWW 




8015:00 


49 FLAG 

DFB 

$00 

{NOT USED 

8016:00 


50 STATUS 

DFB 

$00 

{NOT USED 

8017:00 


51 BUFL 

DFB 

$00 


8018: 


52 BUFH 

DS 

1 


8019: 


53 SAMPLE 

DS 

1 


80 1A: 


54 TMPBLD 

DS 

1 

{USED FOR TEMPORARILY HOLDING A BYTE. 

801B: 


55 ; *»** ft*************** ************ ***************** ******** 


801B: 
8 >HB: 
801B; 
80 IB: 
801B: 
80 IB: 
801B: 


56 

57 

58 

59 

60 
61 
62 


EQUATK 8 

♦**♦**** * »********** «» »»»»*«»*»***********»**»***************** 

THE BASE ADDRESS FOR ALL APPLE SLOTS IS $COOO. 

THE J.B. CARD ADDRESSES BTART AT CSOO WHERE "S" IS THE J.B.'S SLOT NUMBER. 
THE J.B. CARD USES THE ROM LOCATIONS, NOT THE DEDICATED LOCATIONS, ON 
THE APPLE I/O MAP. THE A/D CARD U 8 E 8 THE DEDICATED ADDRESSES SO THE 
A/D BASE ADDRESS IS $C080. 
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0004 : 


63 BELSLT 

EQU 

4 

jJ.B. CARD IN SLOT 4. 

C400: 


64 JBBASE 

EQU 

BELSLT+2 56 +49 1 52 ;THE BASE ADDRESS OF THE J.B. CARD. 

OOJE: 


65 BET HIT 

EQU 

14 

JINTERUPT ENABLE ADDRESS REGISTER ON THE J.B. CARD. 

0009: 


66 CLRIRQ 

EQU 

9 

JJ.B. CARD TIMER2 1NTERUPT CLEAR AND COUNTER START. 

0008: 


67 TMADDR 

EQU 

8 

;LOW BYTE OF TIHBR2 LATCH. 

8010: 


68 CNTL 

EQU 

COUNT 

(LOW BYTE OF ADDRESS CONTAINING THE NUMBER 

801D: 


69 ; 



OF SAMPLES TO BE TAKEN. 

0002: 


70 ODRA 

EQU 

2 

JDATA DIRECTION REGISTER A ON JOHN BELL CARD. 

0003: 


71 DDRB 

EQU 

3 

(DATA DIRECTION RECISTER B ON JOHN BELL CARD. 

8011: 


72 CHTtt 

EQU 

CHTL+1 

jHICH BYTE OF SAME ADDRESS. 

OOOB: 


73 ACR 

EQC 

11 

JAUX REGISTER ADDRESS ON J.B. CARD. 

C408: 


74 T2CL 

EQU 

JBBASE+THADDR jABSOLUTE ADDRESS OF TIMER 2 COUNTER LOW BYTE, 

801 B: 


75 } 



WRITING T2CL WRITES TO THE LATCH. READING T2CL 

801B: 


76 j 



READS FROM THE COUNTER LOW BYTE. 

C409 i 


77 T2CH 

EQU 

JBBASE+CLRIRQ jABSOLUTE ADDRESS OF HIGH COUNTER BYTE. 

801B: 


78 ; 



THERE IS NO TIMER 2 HIGH LATCH BYTE, 

OOOD: 


79 XFR 

EQ8 

13 

J INTERRUPT FLAG REGISTER ADDRESS ON J.B. CARD. 

C40D: 


80 JBIFR 

EQU 

JBBASE+IFR jABSOLUTE ADDRESS OF IFR. 

OOOC: 


81 PCR 

EQU 

12 

J ADORES B OF PERIPHERAL CONTROL REGISTER. 

C40D: 


82 JBFCR 

EQU 

JBBASE+ZFR JABSOLUTE ADDRESS OF J.B* PCR 

0012: 


83 DGPAUS 

EQU 

18 

JTHIS CONTROLS THE DIGITIZING PERIOD. SEE IHNDLR . 

C402: 


84 JBDDRA 

EQU 

JBBASE+DDRA JABBOLUTE ADDRESS OF DATA DIRECTION REGISTER A 

C403: 


85 JilDDRB 

EQU 

JBBASE+DDRB JABSOLUTE ADDRESS OF DATA DIRECTION RECISTER B 

801B: 


86 ; 




801B: 


B7 •++*+****+****«**+*+**+*+«**++*+*+*+**++*+** ******************** 

801B: 


88 i INTERRUPT HANDLER 


801B: 


S9 ; **++*++ ++++*+•+**+*+**+*+**+*+>'*++****+****+*+★**++ A****************** 

801B: 


90 ; 


THIS 

ROUTINE IHieDlATELY RESTARTS SHE J.B. CARD 

801B: 


91 ; 


COUNTDOWN. IT THEN TAKES A SAMPLE FROM THE A/D CARD, 

801 B: 


92 ; 


FIGURES OUT WHERE THE NEXT FREE BUFFER ADDRESS IS , 

801B; 


93 i 


STORES THE SAMPLE THERE, INCREMENTS THE VALUE OF THE 

801B: 


94 i 


NEXT 

FREE BUFFER ADDRESS, AND DECREMENTS THE NUMBER 

801B: 


95 i 


OF SAMPLES YET TO BE TAKEN. 

801B: 


96 j 




801B: 


97 i 




801B:A5 45 


98 IHNULR 

LOA 

$45 

JTHE APPLE INTERRUPT RESPONSE 

8011): 


99 i 



PUTS THE ACCUMULATOR AT $45. 

8010:48 


100 

PHA 



801E: 


101 J 



GET THE ACCUMULATOR FROM $45 

801E: 


102 ; 



AND PUT IT ON THE STACK. 

801E:8A 


103 

TXA 



801F:48 


104 

FHA 



8020:98 


105 

TYA 



8021:48 


106 

PHA 



8022: 


107 ; 




8022: 


108 ; 



RELOAD AND RESET TIHER2 LATCH OH J.3. CARD. 

8022: 


109 j 




8022:18 


110 

CLC 



8023 :AD 08 

C4 

111 

LOA 

T2CL 

JLOAD ACCUMULATOR WITH THE LOW COUNTER BYTE, 

B026 : 


112 ; WHICH BY 

THIS POINT HAS DECREMENTED BELOW ZEF-J. 

8026:69 12 


113 

ADC 

0DGPAU8 

JADD 18 TO IT. SINCE TIHER2 RESTART IS 20 CYCLES 

8028: 


114 J 



AWAY AND THE REFERBICE PUISE5 ARE 10 USEC LONG, 

8028; 


115 S 



TWO REFERBICE PULSES WILL HAVE PASSED 

6028: 


116 ; 



PASSED BY THE TIME TIHER2 IS RESTARTED 

80 28 : 8D 08 

C4 

117 

8TA 

T2CL 

J REMEMBER THAT WRITING T2CL WRITES 

802B: 


118 ; 



TO THE LATCH NOT THE COUNTER. 

802B: 


119 ; 



THIS STEP RELOADS THE TIMER LATCH APPROPRIATELY. 

802B:A9 00 


120 

LOA 

#00 


802D:BA 


121 

NOP 



802E : EA 


122 

NOP 



80?F:EA 


123 

NOP 



8030 :EA 


124 

NOP 



8031:80 09 

C4 

125 

STA 

T2CH 

; ST ORE 0 IN THE HIGH COUNTER BYTE, STARTING TIMER2 

8034 :AE OB 80 

126 

LDX 

CHNL 


8037: BO 80 

CO 

127 

LOA 

$C080 ,X 

JREAD A SAMPLE AND START ON A NEW READING. 

803A:EA 


128 

NOP 



803B:EA 


129 

NOP 



803C :EA 


130 

NOP 



803D;EA 


131 

NOP 



803E:EA 


132 

NOP 



803F ;BD 80 

CO 

133 

LOA 

$C080 ,X 

JREAD SAMPLE STARTED BY PREVIOUS LINE. 

8042 :8D 19 

80 

134 

STA 

SAMPLE 

j STORE THE SAMPLE IN A TEMPORARY LOCATION. 


8045 :AD 17 

80 

135 

LOA 

BUFL 

.LOAD ACCULKULATOR WITH TOE LOW BYTE OF THE 

8048: 


136 ; 



AVAILABLE BUFFER SPACE. 

8048 :8d 55 80 

137 

8TA 

PATCH* 1 

{STORE IT AT ADDRESS P.*.TCH*1. NOTE THAT THIS 

8041)! 


138 i 



ALTERS TOE PROGRAM ITSELF. 

804B:AD 18 

80 

139 

LOA 

BUF1! 

{LOAD THE HIGH BYTE OP THE BOTTOM OF AVAILABLE 

804E:8t) 56 

80 

140 

ST A 

PATCH*2 

{BUFFER SPACE STORE IT AT ADDRESS PATCH* 2. 

8051iAD 19 

80 

141 

LDA 

SAMPLE 

{GET THE STORE SAMPLE OUT OF TOE TEMPORARY LOCATION 

80 54: an 00 

01 

142 PATCH 

8TA 

$0100 

{STORE ACCUMULATOR IN THE ADDRESS GIVEN BY THE 

805/: 


143 ; 



NEXT TWO BYTES. THESE ARE TOE VALUES PUT IN BY 

80 57: 


144 ; 



THE PR7IVIOU6 "STA PATCH" LINES. 

8057: 


145 { 



THE VALUE OP $0100 IS JUST A PLACE HOLDER. 

8057 SEE 17 

80 

146 

INC 

BUFL 

{INCREMENT THE LOU BYTE OF THE ADDRESS OF THE 

80 5A: 


147 i 



LOWEST AVAILABLE BUFFER SPACE TO 

8054! 


148 ; 



GET READY FOR THE NEXT SAMPLE. 

8O5A:D0 03 


149 

BNE 

DECCNT 

{IF BUFL WAS INCREMENTED TO ZERO BUFH HUBT 

80 5C: 


150 ; 



BE INCREMENTED ONE TO TAKE CARE OF THE CARRY. 

80 5C ;EE 18 60 

151 

INC 

BUFH 


805F:CE 10 

80 

152 DECCNT 

DEC 

CNIL 

{DECREMENT THE NUMBER OF SAMPLES YET TO BE TAKEN. 

8062:00 08 


153 

BNE 

LEAVE 

{IF CNTL WAS NOT DECREMENTED TO ZERO 

8064: 


154 { 



THEN LEAVE IHNDLR. 

8064: AD 11 

80 

155 

LDA 

CNTH 

{OTHERWISE CHECK CNTH TO SEE IF IT IS ZERO. 

8067! 


156 ; 



IF CNTH IS ZERO THEN ALL 

8067: 


157 i 



REQUIRED SAMPLES RAVE BEEN TAKEN 

8067 !FO 70 


158 

BEQ 

EXIT 

{IF ALL REQUIRED SAMPLES ARE IN 

8069: 


159 ; 



EXIT FROM THE DIGITIZER. 

8069 :CE 11 

80 

160 

DEC 

CNTH 

{IF CNTH WAS NOT ZERO, DECREMENT IT SINCE 

B06C: 


161 ; 



CNTL WENT THROUGH ZERO. 

8060:68 


162 LEAVE 

FLA 



8060 :A8 


163 

TAY 



806E:68 


164 

FLA 



806F :AA 


165 

TAX 



80/0:68 


166 

FLA 



8071:40 


167 

RTI 



8072: 


168 * *****■>*«*« < «*»*«.■<■»♦**»****-»********♦***********«♦♦**«♦********♦ 

8072: 


169 1 




8072: 


170 ; INITIALIZE 

INTERRUPTS , ETC. 

8072: 


171 ; 




8072: 


172 j ** ********** ************ ********************************** ***** 

8072:48 


173 IN IT 

PHA 



8073 :8A 


174 

TXA 



8074:48 


175 

PHA 



8075:98 


176 

TYA 



8076:48 


177 

PHA 



80 77: AD 11 

80 

178 

LDA 

CNTH 

{MAKE SURE THAT AT LEAST ONE SAMPLE WILL BE TAKEN. 

80 /A: DO 08 


179 

BNE 

OKAY 


807C :AD 10 

80 

180 

LDA 

CNTL 


80/F :D0 03 


181 

BNE 

OKAY 


8081 :EE 10 

80 

182 

INC 

CNTL 

{IF NOT FORCE AT LEAST ONE SAMPLE TO BE TAKEN. 

8084:AD OF 80 

183 OKAY 

LDA 

BUFFER 

{LOAD AND STORE HIGH B'.XE OF DATA BUFFER 

8087 ;8D 18 80 

184 

8TA 

BUFH 

;AT BUFH. 

808A:A9 00 


185 

LDA 

#$00 

{LOW BYTE OF DATA BUFFER IS ASSUMED TO BE ZERO. 

808C:8D 17 

80 

186 

STA 

BUFL 


808F :A9 IB 


187 

LOA 

#IHNDLB 

{PUT LOW BYTE OF IHNDLR ADDR INTO ACC. 

8091 :8D FE 03 

188 

STA 

$3FE 

{STORE AT LOW BYTE OF INTERRUPT VECTOR. 

8094 :AD 14 80 

189 

LDA 

INTADH 

{PUT HIGH BYTE INTO ACCUMULATOR. 

8097: 


190 5 



FOR INFO ON INTADH SEE 

8097: 


191 { 



"DATA STOHACE LOCATIONS" SECTION. 

8097 :8D FF 03 

192 

STA 

$3FF 

{STORE AT THE INTERRUPT VECTOR HIGH BYTE. 

809A: 


193 J 




809A: 


194 ; 




80 9A: 


195 ; 




809A: 


196 -.SELECT 

TIMER2 FUNCTION ON J.B. CARO BY SETTING UP THE AUX. CTRL. REG. 

809A: 


197 ;SET UP TOE 

INTERUPT ENABLE REGISTER TO ALLOW ONLY TIMER2 INTERRUPTS. 

80'iA;A2 OB 


198 

LDX 

tfACR 

{J.B. ADDRESS OF AUXILIARY CONTROL REGISTER. 

8C9C, -\9 20 


199 

LOA 

#$20 

{TO BE WRITTEN TO TOE AUX CONTROL REGISTER. 

80311:3.9 00 C4 

200 

STA 

JBBASE.X 

{SELECTS TIMER2 EXTERMAL COUNT MODE. 

80Al:Afi OK 


201 

LDX 

#SETINT 

;J.S, ADDRESS OF THE INTERRUPT ENABLE REGISTER. 

80A3 :A9 AO 


202 

LOA 

#$A0 

{CODE TO ENABLE TIMER2 INTERRUPT. 

80 A5 : 90 00 

C4 

203 

!5TA 

JBBASE.X 

'{STORE CODE IN INTERRUPT ENABLE REGISTER. 

80A8 :A9 5F 


204 

LDA 

#$5P 

{CODE WHICH DISABLES ALL OTHER INTERRUPTS. 

80AA:9D 00 C4 

205 

STA 

JBBASE.X 

{DISABLE ALL OTHER J.B. INTERRUPTS • 


80 AD :A9 
80AF:8P 
80 B2: 

80 B2: 

80 B2 s 
80B2 : 
80B2 : 
80B2 ; 

80 B2: 

80 B2; 
80B2 : 

80 B2; 

80B2-.A2 

80B4:A9 

80B6-.9D 

80B9:A2 

80BB:A9 

8(JBD :9D 

80C0: 

80C0:AD 

80C3:P0 

80C5: 

80C5: 

80C5: 

80C5 :A2 

8007 !A9 

80C9:9D 

80CC: 

80CC: 

80CC : 

80CC: 

80CC :A2 

BOCEiAB 

80D0:9D 

80D3;A2 

80D5;A9 

80D7:9D 

SODA : 

80DA: 

80 DA: 

SODA: 

SODA: 

80DA:68 

80DB :A8 

80DC:66 

80DD:AA 

80DE:66 

80DF:S8 

80EO:6O 

80E1: 

80 El: 
80E1: 

80 El: 

80E1:48 

80E2 :8A 

80E3-.48 

80E4:9» 

80E5-.48 

80 E6: 

80 E6: 
80E6: 
80E6: 
80E6: 
80E6 : 
80E6: 
80E6 ;A9 
60E8 :A2 
S0EA:9D 
80ED :A9 
80EF; 
80EF:8D 


00 

03 C4 


206 

207 

206 

209 

210 
211 
212 

213 

214 

215 


LDA 

STA 


800 

J8DDRB 


{CONFIGURE DATA I/O P.ECIBTEB B LIMES A3 INPUTS. 


FINISH UP INITIALIZATION UPON RECEIVING A POSITIVE TRANSITION ON LIME CA1 
FIRST CLEAR THE IFF. THEM 8ET UP THE PERIPHERAL CONTROL REGISTER TO 
,SET THE CA1 INTERRUPT FLAG OH A NEGATIVE CA1 TRANSITION. THEM SIT IN A 
WAIT LOOP, MONITORING THE IFR.ONCE A FLAG 18 SET START THE COUNTER AND 
START DIGITIZING. 

NOTE THAT THE CAl INTERRUPT ENABLE BIT IS t INOTt I SET. 



216 

-,1T IS 

ONLY NECESSARY TO DETECT IPR BEING SET. 


217 

;D0 NOT WANT 

INTERRUPT 

SEQUENCE TO START. 

OD 

218 


LDX 

Fife 

FF 

219 


LDA 

#$PF 


00 C4 

220 


STA 

JBBASS.X 

{CLEAR THE IFR. 

OC 

221 


LDX 

#PCR 


00 

222 


LDA, 

#00 

{INTERRUPT SET ON A NEGATIVE TRANSITION ON CAl . 

00 C4 

223 


STA 

JBRASE.X 

{WRITE TO PCR. 


224 

S 




OD C4 

225 

BTRTDG 

LDA 

JBIFR 


FB 

226 


BEQ 

STRTDG 

{IF IFR IS CLEAR, STARTING SIGNAL HAS 


227 

> 



NOT YET BEQ! RECEIVED. 


228 

» 



IF IFR NOT CLEAR THEN CONTINUE TO AWAIT TRIGGER. 


229 

* 




OD 

230 


LDX 

#IPR 

{CLEAR THE IFR 

FP 

231 


LDA 

#8FF 


00 C4 

232 


STA 

JBBASE.X 



233 

* 





234 

a LOAD UP TIMER2 LATCH 



235 

♦ 

9 





236 

f 




08 

237 


LDX 

FTMADDR 


12 

238 


LDA 

DGPAUS 


00 C4 

239 


STA 

JBBASE.X 


09 

240 


LEX 

Fclrirq 


00 

241 


LDA 

#800 


00 C4 

242 


STA 

JBBASE.X 

{START COUNTER OH TIMER2. 


243 

• 




244 

J RESTORE PEG1STEKS 



245 

{SINCE 

THE J. 

,B. COUNTER 

IS RUNNING AT THIS POINT, THE FOLLOWING STEPS 


246 

{SHOULD NOT TAKE LONGER 

THAN THE TOTAL TIME COUNTED DOWN BY THE J.B. CARI 


247 

t 





248 


PLA 




249 


TAY 




250 


PLA 




251 


TAX 




252 


PLA 




253 


CLI 




254 

255 

:**+**i 

RTS 

I********** 



256 

257 

258 ; 

259 FINISH 


FINISH UP AND DISABLE J.B. CARD. 


260 

261 

262 

263 

264 

265 

266 

267 

268 

269 

270 


{SAVE THE REGISTERS 


PHA 
TXA 
PHA 
TYA 
PHA 

CAUTION! (Ill 

THE FOLLOWING ROUTINE IS REACHED BY A JUMP PROM IHNDLR. EXIT'S PURPOSE 
IS TO RETURN TO THE BASIC PROGRAM FROM THE GALL TO THE WAIT LOOP. 

TO DO THIS THE ADDRES8 OF THE 

WAIT LOOP, PLACED ON THE TOP OF THE STACK BY THE INTERRUPT, 

MUST BE SKIPPED OVER, BRINGING UP THE ADDRESS OF THE BASIC STATEMENT 
AFTER THE CALL TO THE WAIT LOOP. 


7P 

271 EXIT 

LDA 

#S7F 


OE 

272 

LDX 

#SETINT 

{DISABLE J.B. INTERRUPTS. 

00 C4 

273 

STA 

JBBASE.X 


65 

274 

LDA 

#865 

{THESE FOUR STEPS RESTORE 


275 { 



ORIGINAL INTERRUPT VECTOR 

FF, 03 

276 

STA 

83 FE 
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80?2:A9 PP 

277 


LDA 

#8ff 


80F4:8D PP 03 

278 


8TA 

83PF 


80P7I68 

279 


PU 


{RESTORE TOE 

B0P8:AB 

280 


TAV 



80P9 :68 

281 


PU 


{REGISTERS. 

80FASAA 

282 


TAX 



80 FB;68 

283 


PU 



80 PC :8D 1A 80 

284 


BTA 

TMPHLD 

{SAVE ACCUHDUTOR TEMPORARILY 

80FF:28 

285 


PLP 


{TAKE PROCESSOR STATUS BYTE PROM STACK. 

8100: 

286 

» 



THIS MIMICK8 

8100; 

287 




TOE RTI INSTRUCTION WITHOUT RETURNING TO TOE 

8100; 

288 

* 

| 



INTERRUPT CALL 

8100:68 

289 


PU 


{WASTE TOE RETURN- TD-WA IT” PROM- INTERRUPT ADDRESSES 

8101:68 

290 


PLA 



8102:AD 1A 80 

291 


LDA 

TXP11LD 

{GET ACCUKUUTOR FROM TEMPORARY HOLDING SPOT. 

8105:S8 

292 


CL I 



8106:60 

293 


RTS 


{RETURN TO BASIC PROGRAM, NOT TO WAIT LOOP 

8107:58 

294 WAIT 

CL1 



8108 :EA 

295 

BACK 

HOP 



8109:40 08 81 

296 


JMP 

BACK 


810C: 

297 

l 




8100: 

298 

t 




8100: 

299 

* 
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II .2 TEXTIFY 


40 ROC PROCRAM TEXTIFY. RODNEY BTOLTZFUS, PEB.l 1984. 

50 RBI THIS PROGRAM CONVERTS A BINARY FILE TO A TEXT FILE 

51 REM WITH THE FORMAT XYZXYZXYZ WHERE XYZ IB A THREE PLACE DECIMAL 

INTEGER. 

52 REM THE USER SUPPLIES THE STARTHC ADDRESS OF THE BINAY FILE IN 

APPLE MEMORY. 

53 REM THE RECORDS OF THE C OPERATED TEXT FILE HAVE 20 THREE DIGIT 

NUMBERS . 

100 D$ - CHRS (4) 

110 INPUT "AT WHAT ADDRESS DOES THE BINARY DATA START IN MEMORY? START 
120 INPUT "HOW MANY BYTES TO BE CONVERTED? SHOULD BE A MULTIPLE OF 20."; 
LNG 

130 INPUT "NUMBER OF FILE (#) ? (FILE NAME IS ROCKET TKXT#)";K$ 

190 PRINT "RUNNING..." 

200 PRINT D$”OPEN ROCKET TEXT";N$ 

300 PRINT D$"WRITE ROCKET TEXT";N$ 

410 RBI START GIVES THE ADDRESS OF THE FIRST BYTE TO BE TRANSMITTED. 

500 FOR J - 1 TO (INC / 20) 

550 C$ - "" 

600 FOR I « 1 TO 20 

610 RBI ADD 1000 TO EACH NUMBER TO GIVE EACH NUMBER THE SAME LENGTH. 
700 X » 1000 ♦ PEEK (START + (J - 1) * 20 ♦ (I - 1)) 

750 REM BUILD STRING WITH RIGHT 3 CHARACTERS OF THE OFFSET DATA 
800 C$ - C$ ♦ RIGHTS (( SIRS (X)),3) 

900 NEXT I 
910 PRINT C$ 

920 Y - FRE (0) 

1000 NEXT J 

1100 PRINT D$"CLOSB ROCKET TEXT";N$ 

1110 PRINT "FINISHED! DATA ON DISK" 

1200 END 


II. 3 APL2CDC 


PROGRAM APL2CDC(INPUT, OUTPUT, APLFIl.,AP).OUT,TAPEl -INPUT, 

* TAPE 2 “OUTPUT, TAPE3-APLFIL,TAPE4-AP LOUT) 

C 

0 THIS PROGRAM CONVERTS AN ICE PILE TO A FILE CAPABLE OF 

C BEING READ BY THE PROGRAM SPECTRA AND OTHERS WITH THE 

C SAME FORMAT NEED. THE ICE FILE CONVERTED IS TRANSFERRED 

C FROM THE »J>PLE (SEE APPLE PROGRAM TEXTIFY FOR FORMAT). 

C THE APPLE FILE IS CALLED ROCKET TEXT. 

C THIS PROGRAM ASSUMES A 5KHZ DIGITIZING RATE AND A .05 SECOND 
C RECORD LENGTH. TO CHANGE THIS, CHANGE THE LENGTH OF THE ARRAYS. 

C AND FIX THE TIME SPACING BETWEEN RECORDS (I.E. CHANCE RECPRD, 

C THE TIME PERIOD OF ONE RECORD) 

C THE OUTPUT RECORDS HAVE THE FORMAT [F/DATA/F/F/F] WHERE F IS 
C FILLER (I.E. IS ONLY TAKING UP SPACE) AND DATA IS THE DATUM. 

C THIS FIVE POINT SEQUENCE IS CALLED A FRAME AND THERE ARE 250 

C FRAMES PER .05 SECOND RECORD. 

0 TEXTIFY AND APL2CDC WERE CREATED BY RODNEY STOLTZFUS 

C DECEMBER, 1983 

C AS OF FEB. 3, 1984 THE FORMAT OF THE APPLE FILE IS 20 

C 3 DIGIT NUMBERS PER TEXT LINE. THIS CAUSES PROBLEMS AS THE 

C APL2CCC OUTPUT SHOULD BE IN BLOCKS OF 250. THE SOLUTION IS 

C TO READ IN 500 POINTS AT A TIME AND PRINT OUT 2 ARRAYS AT 

C ONE PASS THROUGH SUBROUTINE REFORM. ESSENTIALLY THIS MEANS 

C THAT THE NUMBER OF INPUT POINTS SHOULD BE A MULTIPLE OF 500 

C RATHER THAN 250. 

C 

C 

c********* ******** *********************************************** 

c 

C VARIABLE DICTIONARY 

C 
C 

C I ARRAY: 

C 
C 

C ARRAY : 

C 
C 

C DATNUM: 

C 

C DATA : 

C 

c I, j,k : 

c 

C BLOCKS : 

c 

S RECPRD: 

c 

C*************** *************************************************** 

c 

INTEGER DATA ( 500 ) , I , J , K , DATNUM, BLOCKS , J1 9 
REAL TIME 
C 
C 

PRINT*, "WHAT IS THE UNIVERSAL TIME OF THE FIRST RECORD?" 
READ*,TIME 

PRINT*, "HOW MANY DATA POINTS ARE IN THE INPUT FILE?" 

PRINT*, "(MUST BE MULTIPLE OF 500. SEE OPENING PROGRAM COMMENTS)" 
READ* .DATNUH 
C 

C READ IN THE DATA IN DOUBLE BLOCKS 
C 

BLOCKS-DATNUM/250 
DO 110 1-1, BLOCKS, 2 
DO 100 J-l, 500,20 
J19»J*19 

READ (3,900) (DATA(K) ,K“J, J19) 

100 CONTINUE 
C 

C INSERT THE FILLER DATA AND THE TIME CODE FOR TWO BLOCKS. 


HOLDS THE INTEGER DATA BEFORE TRANSMITTING. I ARRAY 
IS EQUIVALENT TO ARRAY (2 TO 1251). 

HOLDS THE 1 ARRAY AND THE TIME CODE. THE TIME CODE IS 
IN ARRAY ( 1 ) . 

IS THE TOTAL NUHBER OF DATA POINTS IN THE FILE 
IS AN ARRAY THAT HOLDS 500 PIECES OF RAW DATA. 

GENERAL PURPOSE COUNTERS. 

GIVES THE NUHBER OF 250 POINT BLOCKS IN THE INPUT FILE. 
GIVES THE TIHE DURATION OF ONE RECORD. 


n n o s n r n o o o f) o o 


WRITE THEM OUT TO TAPE4 (8EE PROGRAM STATEMENT) , 

CALL REF0RH(DATA,I,TIHE) 

110 .CONTINUE 
900 FORMAT (20 13) 

STOP 
END 

SUBROUTINE REFORM INSERTS THE EXTRA DATA POINTS, PUTS ON 
THE TIME MARK, PUTS THE DATA INTO ARRAY, AND WRITES OUT ARRAY. 
IT PUTS OUT TWO ARRAYS IN A PASS. 

SUBROUTINE REFORM(DATA,RECNUH,TIME ) 

REAL ARRAY (1251), TIME, RECPRD 
INTEGER RECHUH, DATA( 500 ) , IARRAYU 250 ) , K 
EQUIVALENCE (ARRAY(2), lARRAY(l)) 

PUT IN EXTRA GARBAGE DATA POINT8 


SET THE RECORD DURATION TIKE. FOR 5 KHZ IT SHOULD BE .05 SECONDS. 

RECPRD-. 05 
DO 210 K-1,2 
DO 200 J-1,250 
I ARRAY ( (J-l)*5+l) “9999 
I ARRAY ( ((J-l )*5)+2)-DATA(J+250*(K-l) ) 

IARRAY({(J-l)*5)*3)-9999 
IARRAY( (( J-l)*5)+4)-9999 
IARHAY( ((J-l)*5)+5)-9999 
200 CONTINUE 
C 

C THE NEXT STATEMENT ATTACHES A TIME CODE RELATIVE TO TIME. 

C WITH A RECPRD SECOND RECORD LENGTH , BACH AFFIXED TIME CODE SHOULD 
C BE RECPRD SECONDS GREATER THAN THE PRIVIOUS ONE. 

C (RECNUM-1) IS USED SINCE THE NTH RECORD SHOULD HAVE A TIME CODE 
C OF (N-l)*RECPRD. THE (R-l) IS USED SINCE SUBROUTINE REFORM PUTS 
C OUT 2 RECORDS, ONE AT K-l AND ONE AT K-2. 

C 

ARRAY U ) “TIME* (FLOAT( (RECNUM-1 )+(K-l ) )*RECPRD) 

WRITE (4) ARRAY 
210 CONTINUE 
RETURN 
END 
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PROGRAM SPECTRA IS All INTER- ACTIVE ROUTINE WHICH ALLOWS VUE 
COMPUTATION OF POWER SFECTRA TO DE USED IN THE CALCULATION 
OP FINE STRUCTURE IRREGULARITIES . THE FINE STRUCTURE DATA 
FILE TO BE ANALYZED IS ENTERED THROUGH DEVICE 3 . THE USER 
IS THEN PROMPTED TO PROVIDE INFORMATION REGARDING TOE DATA 
AND THE ANALYSIS. A LINEAR FIT MAY ALSO BE PERFORMED ON 
TOE DATA AND TOE SPECTRAL INDEX CALCULATED, THE USER HAY 
STORE VALUES OF THE SPECTRAL INDEX AND MAGNITUDE IN ARRAY 
INDEX. PROGRAM 8PINDEX HAY THEN BE USED TO PLOT THE 
SPECTRAL INDEX AND MAGNITUDE AS A FUNCTION OF ALTITUDE. 

PRIOR TO EXECUTION OF SPECTRA, TOE USER MUST PROVIDE 
CERTAIN PARAMETERS. THESE PARAMETERS ARE 

1) CAIN — THIS IS TOE CAIN OF THE EXPERIMENT. 

2) CHANNEL — THIS REFERS TO TOE CHANNEL NUMBER 

OF TOE REQUIRED DATA AS STORED ON 
A CDC FORMAT TAPE. 

3) MIN HAG — THIS IS THE LOWER BOUNDARY ON THE 

PLOTTING AXIS. 

4) BAHDCAL — THIS IS TON ARRAY WHICH CONTAINS 

THE DATA CALIBRATIONS AS DETERMINED 
PR CM PROGRAM CWALCAL. 

5) LOLIMIT — LOLIHIT IS THE LOWER BOUND OH THE 

FREQUENCY OF DATA POINTS TO BE INCLUDED 
IN THE LEAST SQUARES FIT. 

6) HILIMIT — HILIMIT IS THE UPPER BOUND ON TOE 

FREQUENCY OF DATA POINTS TO BE INCLUDED 
IN TOE LEAST SQUARES FIT, 

7) FC -- THIS IS THE UPPER 3DB FREQUENCY OF THE 

TELEMETRY LOW PASS FILTER. 

8) FAC — THIS IS THE UPPER 3dB FREQUHICY OF THE 

AC AMPLIFIER. 

9) FLE — THIS IS THE UPPER 3D3 FREQUENCY OF THE 

LOGARITHMIC ILEC7R0METER. 

10) BBAND — THIS IS THE HUMBER OF THE BROADBAND 

AMPLIFIER OUTPUT BEINC PROCESSED.* (TO 
BYPASS DEBIPH, SET BBAND .EQ, 0) 

THESE VALUES ARE QJTERED'UNDER TOE TITLE "EXPERIMENT 
PARAMETERS" IN THE PROGRAM. (WRITTEN BY BRUCE TOMEI, 1983). 

THIS VERSION OF SPECTRA IS MODIFIED FROM THE ORIGINAL VERSION 
BY BRUCE TOMEI. AN ATTEMPT HAS BEEN MADE TO DOCUMENT 
MAJOR CHANGES. DEBUGGING OF TOE ORIGINAL 18 NOT 
RECORDED. THERE IS AN EXTENDED 

PLOTTER ROUTINE ALLOWING INCLUSION OF THE ROCKET SHOT 
NUMBER, THE LAUNCH TIME (IN UNIVERSAL TIME), THE NUHBER 
OP SECONDS INTO TOE LAUNCH AT WHICH SPECTRA ARE DESIRED, 

THE ROCKFT SPEED AT THAT TIME, THE ALTITUDE, THE ROCKET 
ANGLE OFF VERTICAL, THE FFT BLOCK LENCTO, THE SEGHENT 
OVERLAP, THE DIGITIZATION RATE, THE ROCKET SPIN RATE, 

AND THE CAIN OF THE EXPERIMENT. 

NEW VARIABLES AND CONSTANTS ARE: SHQTNUM,ALTITUD, AND ROCK6PD. 

THIS VERSION OF SPECTRA DOES NOT ASK FOR TOE ANGLE, 

OR THE SPIN RATE OF THE ROCKET. 

RODNEY STOLTZFUS. 10/11/83 

ANOTHER MODIFICATION: 

SPECTRA USED TO PLOT A STRING OF SEQUENTIAL PLOTS. YOU COULD 
PLOT AS MANY AS YOU WISHED OUT OF THE AVAILABLE SEGMENTS . 
HOWEVER, YOU COULD NOT PICK OUT INDIVIDUAL PLOTS. FOR EXAMPLE, 
YOU COULD NOT PLOT ONLY THE FIFTH SEGHENT; YOU HAD TO 
PLOT SEGMENTS 1,2, 3, 4, AND 5. 

THE MODIFICATION ALLOWS YOU TO CHOOSE ANY AVAILABLE PLOT OR 
PLOTS . THE COMPUTED AVERAGE IS OVER THESE SELECTED 
PLOTS RATHER THAN OVER THE WHOLE RANGE OF SEGMENTS. 

THIS ADMITTEDLY CONSTRAINS THE AVERAGE PLOT TO COVER ONLY 
THOSE SPECTRA WHICH ARE ACTUALLY PRINTED. YOU MUST PLOT 
ALL SPECTRA YOU WISH TO HAVE INCLUDED IN THE AVERAGE. 

THIS MODIFICATION IS PERFORMED BY PUTTING TOE DESIRED PLOT NUMBERS 
IN AN ARRAY CALLED PLOTS . YOU CAN INCLUDE UP TO 99 PLOTS . 
SUBROUTINE SELECT HANDLES TOE CONSTRUCTION OF TOE ARRAY CALLED 
PLOTB. 

2/16/84 RODNEY STOLTZFUS 


126 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

G 

c 


c 

c 

c 


c 

c 

c 


MODIFICATION: SPECTRA HOW HAS THE ABILITY TO SAVE AVERAGE 

PLOTS FOR UTER PROCESSING, YOU CAN PUT AN AVERAGE PLOT 
INTO MEMORY LOCATIONS (CALLED PLOT BAYS) AND CO BACK TO 
CREATE MORE AVERAGE PLOTS . 

PROCESSING INCLUDES AVERAGING TWO PLOTS, TAKING TOE DIPPER BICE 
OP TWO PLOTS , PERFORMING AN INCREMENTED WINDOW AVERAGE, 
AND PLOTTING THE FINAL RESULT. 

THERE ARE CURRENTLY TWO PLOT BAYS PLUS A WORKING BAY IN WHICH 
PROCESSED DATA IS KEPT. IP NECESSARY, MORE BAYS CAN BE 
IMPLEMENTED BY DECLARING NEW ARRAYS AND ADJUSTING TOE 
SUBROUTINE CALLS. 

A PLOT BAY CONTAINS THREE ARRAYS: fHE FREQUENCIES Op THE PPT, 
THE CORRESPONDING MAGNITUDES, AND ASSOCIATED PARAMETERS. 
THE PARAMETERS ARE THE TUB OF THE PLOT, THE SIZE OF TOE 
DATA ARRAYS, AND THE FFT SIZE. 

NOTE: DEPENDING ON TOE NUMBER OF PLOT BAYS, THE PROGRAM MAY 
REQUIRE MORE CENTRAL MEMORY THAN IS USUALLY GIVEN TO 
A RESEARCH ACCOUNT. IF THE PROGRAM WILL NOT RUN AT 
NORMAL TIMES , YOU MAY NF.ED TO HAVE YOUR CENTRAL MEMORY 
LIMITS EXPANDED. 

3/ 12/84 


SPECTRA WAS WRITTEN BY BRUCE TOMEI, 1983. 

SPECTRA REQUIRES TOE FOLLOWING EXTERNAL INFORMATION: 
THE IMSL LIBRARY (ACCESSED BY GRAB, IKSL) 

THE GCS LIBRARY (ACCESSED BY GRAB,GC8(TEKT,3D) ) 

AN INPUT FILE IN STANDARD ROCKET TILE FORMAT 


PROGRAM SPECTRA (INPUT, OUTPUT, DATAll .INDEX, TAPEl-INPUT, 

* TAPE2 “OUTPUT, TAPE3**DATA1 1 ,TAPE4-INDEX) 

INTEGER N , K , M , OVLP , WINDONO , DIS PLAY , OPT ION , CHANNE L , IBRS , IW1N , S , 

* NUMBER ,FIT1,PIT2, REPEAT, M2, BRAND, DIGHATE, JJ,PLOT8(100) , 

* TOTAL , QUE3 , PLTS IZE , FFTS IZE 

REAL DATA(4096) ,MAG(2049) ,MACAVE(2049 ) ,PR2Q(2049> ,X(4096) , 

* BANDCAL(5) ,MAXKAG,HAXFREQ,LTIME ,STIME,SHPTIME,8ECS .GAIN,! 

* MINHAG.KREAL, SLOPE ,YINT,HILIMIT,LOLIMIT,FC, FAC, FLE.SHOTOUH, 

* ALT ITUD , ROCKS PD , SAVHAG1 ( 102 S > , SAVMAG2 (1025), 

* WORKHAG ( 1 025 ) , S AVFRQ1 ( 1 02 5 ) , SA VFRQ2 ( 1 02 5 ) , 

* W0RKFRQ( 1 025 ) , PARAMSl ( 10 ) , PARAHS2 (10) , 

* WKPARMB(10),AVPAKHB(10) 

COMPLEX Z(2049) 

EXTERNAL CABS 

DATA Z , DATA , MAG , MAG A VE , X, FREQ , SA VHAC1 .8AVMAC2 , WORKHAG, 

* : SAVFRQ1 .SAVFRQ2 .WORKFRQ, PARAMSl , PAR AMS 2 , 

* WKPARMS , AVPARMS / 2.049* ( (0 .0 ,0.0)) ,4096*(0 .0) , 

* 2049*(0 .0) ,2049*(0.0) ,4096*(0 .0) ,2049*(0.0) , 

* 1025*(0 .0) ,1025*(0.0 ) ,1025*(0 .0) , 

* 1025*(0.0) ,1025*(0.0) , 1025* (0.0 ) ,10*(0 .0) , 

* 10*(0.0),10*<0.0),i0*(0.0)/ 

I 

REMINDER TO CHECK EXPERIMENT PARAMETERS. 

PRINT*, "RB1INDRR:MAKE SURE YOU HAVE SET TOE EXPERIMENT" 

PRINT*, "PARAMETERS CORRECTLY (I.E. GAIN .BANDCAL, ETC. )" 

PRINT*," " 

ENTER USER-SPECIFIED VARIABLES 


PRINT*, "ENTER THE LAUNCH TIME. (HH MM SEC.)" 

READ*, IHRS.IMIN, SECS 
LTIME«FLOAT(IHRB*3600+1MIN*60)+SEC3 
10 PRINT*, "ENTER TOE DESIRED STARTING TIME, (SECONDS AFTER LAUNCH)'* 
READ*, ST IMP. 

IF(REPEAT .EQ. l)GOTO 50 

PRINT*, "ENTER TOE NUMBER OF DATA POINTS." 
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20 
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100 

40 


READ*,N 

PRINT*, "ENTER THE BE CHEN T LENGTH (MUST BE AK EVDI INTEGER)" 
REAT*,L 

PRINT*," NUMBER OP SEGMENTS OVERLAP" 

B-N/L 
PRINT4 ,8 ,0 
IPO .EQ, 1) GOTO 40 
00 100 1 - 2 , K 

KREAL-FLOAT<N-<I-t))/FLOAT(Wl-l)) 
ip(kreal-aikt(kreal) 1100,20 

K-IFIX(KHEAL) 

PRINT4,K,I-1 
FORMAT ( 1 3X , 14 , 1 9X , 14 ) 

CONTINUE 

PRINT*, "ENTER HIE NUHBER OP 8ECHENTS " 

PRINT*, "THE DATA SHOULD BE DIVIDED INTO" 

READ*,K 

OVLP-(L*K-N)/(K-l) 

PRINT*, 'ENTER THE DESIRED WINDOW. 

PRINT*," 

PRINT*, " 

PRINT*," 

PRINT*," 

PRINT*," 

READ*, WIN DONO 

PRINT*,"QITER THE PFT BLOCK LENGTH 
READ*,H 
M2-M/2+1 

PRINT*, "ESTER THE DESIRED plOT, 

PRINT*," 

PRINT*," 

READ*, OPT ION 
FIT2-2 
DISPLAY-2 
PIT1-1 

PRINT*, "DO YOU WISH TO PLOT INDIVIDUAL PERIODOGRAMS 7 1-YES" 
PRINT*," 2 -NO" 

READ*, DISPLAY 


1 - RECTANGULAR" 

2- BARTLETT" 

3 - HANNING" 

4- UAMNINO" 

5 - BLACKMAN" 

6- BLACKMAN-UARRI8" 

(MUST BE A POWER OF 2)" 


1- LINEAR" 

2- LOGAR ITHMIC" 

3 - NO PLOT" 


THE FOLLOWING ALLOWS THE USER TO SELECT INDIVIDUAL PLOTS TO 
BE PLOTTED (AND INCLUDED IN THE AVERACE). PLOTS (100) 

HOLDS THE RESPONSE TO THE QUESTION. 

MODIFICATION TO ORIGINAL SPECTRA PROGRAM. SEE INITIAL PAOGBAH 
COMMENTS OR COMMENTS FOR SUBROUTINE SELECT. 2/16/84 


PRINT*, "DO YOU WISH TO DISPLAY AND AVERAGE ONLY CERTAIN PLOTS T" 
PRINT*, ' 1-YES, I ONLY WANT SOME OP THE 8ECHENTB" 

PRINT*," 2 -NO, I WANT ALL SEGMENTS AVERAGED" 

READ*. PLOTS (100) 

C 

C CALL SUBROUTINE SELECT. THIS WILL PUT THE D'tSIRED PLOT NUMBERS 

C INTO THE ARRAY CALLED PLOW . YOU CAN INCLUDE UP TO 99 PLOIB . 

c 

CALL SELECT (PLOTS, NUHBER, K, TOTAL) 

PRINT*, "FERFOBH LINEAR PIT TO INDIVIDUAL PERIODOGRAMS? 1— YES" 
PRINT*," 2—NO" 

READ*, PI T2 

30 PRINT*, "PERFORM LINEAR FIT TO AVERAGE!) PERIODOCRAM? 1— YES" 

PRINT*," 2 — NO " 

READ*, FIJI 

PRINT*, "WHAT 18 THE NUHBER OF THE ROCKET BHOT?" 

RHAB*,BHOTHUH 

C PRIHT*,"WHAT IS THE ALTITUDE FOR THIS ANALYSIS?" 

C READ* , ALTITUD 

C PRINT*, "WHAT IB THE ROCKET SPEED AT THIS ALTITUDE?" 

C READ*,R0CK8FD 

PRINT*, "WHAT IS THE DIGITIZING RATE OF THE DATA?" 

READ*, DIGRATE 

70 PRINT*, "ANY CHANGES? 1-YES" 

PRINT*," 2-NO" 

read*,s 

XF(S .EQ. 1 )COTO 10 


C 
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experiment parameters 

CAIN-35.000 

CHANNEL-2 

HINHAG-O .000000000001 
BAN DCAL(l)— 115.50 
BANDCAL(2)-5 .0600 
BANDCAL(3 >-127.000 
BANDCAL(4)-246.000 
BANDCALC 5 >—370 .000 

PRINT*, "DO YOU WANT TO CHANGE THE BAND CALS? 1-YB8" 

PRINT*," 2-NO" 

READ*,QUES 

IF (QUES .EQ. 2) GOTO 75 

PRINT*, "ENTER THE BANDCALS AS YOU ARE ASKED:" 

PRINT*, "BANDCAK 1)7" 

read*,bandcal(1) 

PRINT*, "BANDCAL<2)7" 

READ*. BANDCAL{ 2) 

PRINT*, "BANDCAL(3)7" 

READ*,BANDCAL(3) 

PRINT*,"BANDCAL(4)7" 

READ*,BANDCAL(.4) 

PRINT*, "DANDCAL(5>?" 

read*,ban;ocal(5) 

'5 LOLIMIT-IOO.O 
HILIMIT-2500 .0 
FC-1050 .0 
FAC-1880,0 
FLE-2500,0 

DBAND ■ 1 IF YOU WANT TO SELECT YOUR OW) LOWER 3DB DEEMPHASIS 
POINT. 

BEAND-2 

BEGIN PROCESSING 
50 REWIND 3 

CALL LOAD (DATA, N , CHANNEL, LTI HE, ST IHE.SHPTIKE) 

CALL CALID(DATA,N, RAN DCAL, GAIN) 

CALL SPECTRH(DATA,KAC,MAGAVE,PREQ,X,Z,N,K,M,M2,L,OVLP,WINIK)NO, 

* SMPT IME , OPT ION , DI8 PLAY , NUMBER ,M INMAG , FI 12 , FIT1 , 

* LOLIHIT ,HILIMIT,FC ,FAC ,FLE,BBAND,LTIME ,8TIME,SHOTNUH 

* DIGRATE, TOTAL, PLOTS) 

CALL FILTGjV?.(MAGAVE,FREQ,M2 ,FC , PAC , PLE > 

CALL DEftt?H(MACAVE ,PREQ ,M2 , BRAND) 

CALL MAXH IN (MAG AVE , FR EQ , MAXMAG ,MAXFREQ,M,M2) 

JJ-0 

SETTING JJ-0 WILL CAUSE 0 TO BE PRINTED ON TOE AVERAGE 
PLOT. 

CALL PLOTTER(MACAVE, FREQ ,M ,M2 , MAXMAG, KIN HAG, MAXFREQ, OPTION , 

* SLOPE , VI NT , FITl , LOL IMIT, H1LI MIT , ST IKE , S HOTNUM, 

* ALTITUD.ROCKSPD.IHRS ,IHIN ,BECS .DIGRATE ,OVLP,CAIN ,JJ) 

PRINT*, "DO YOU WISH TO STORE THE SPECTHAL INDEX? 1— YES" 

PRINT*," 2— NO" 

READ*, 8 

IF(S ,EQ. 2)COTG 60 

CALL STORE ( SLOPE , YINT , LTIKE , ST IME , SMPT 1ME , N ) 


60 PRINT*, "ARE YOU FINISHED? 1-YES" 

PRINT*," 2-NO" 

READ*,S 

IF(S .EQ. l)STOP 

C THE FOLLOWING SECTION ALLOWS FURTHER ANALYSIS OP SPECTRA. 

C SPECIFICALLY IT ALLOWS THE USER TO SAVE THE AVERAGE PLOT MADE BY 
C SUBROUTINE BPECTRM, AND ALLOWS TOE USER TO CREATE MORE SPECTRA. 

C TOE SPECTRA ARE SAVED IN ARRAYS SAVMAG1 AND SAVMAG2 . 

C EACH ARRAY GAN HOLD UP TO 1025 POINTS (CORRESPONDING TO A 2048 
C POINT FFT. 
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C FREQUENCIES ASSOCIATED VITO TOE SAVED HAGNITIUDE6 ARE IN ARRAYS 
C SAVFRQl AND SAVFRQ2. CERTAIN PLOT PARAMETERS ARE STORED IN 
C PARAMS1 AND PARAM82. THE ORDER OF BTORACE 7A TIME, NUMBER OF POINTS, 

C AND FFT SIZE. 

C IN ADDITION TO THESE THREE ARRAYS THERE IS A WORKING PLOT CALLED 
C VORKHAC. WORKFRQ AND WKPARHS STORE TOE FREQUENCIES AND THE 
C ASSOCIATED FREQUENCIES AND PARAMETERS. 

C 

C FIRST, BTORE TOE AVERAGE PLOT IN ONE OF THE PLOT BAYS. ALSO SET UP 
C AN ARRAY AVPARAM . 

C 

PRINT*, "DO YOU Visa TO STORE THE LAST AVERAGE PLOT? 1-YE6" 

PRINT*," 2-NG" 

READ* , QUE8 

IF (QUES .EQ. 2) GOTO 100 
AVPARMS(1)-8TIME 
AVPAP.MS ( 2 ) *M2 
AVFARItS (3) *H 

DISPLAY WHICH PLOTS ARE WHERE 


CALL PRTS TAT ( PAR AKS I , PAR AMS 2 , WKPARHS ) 

PRINT*, "IN WHICU BAY DO YOU WISH TO STORE TOE AVERAGE PLOT" 
PRINT*," 1-BAY#1” 

PRINT*," 2-DAY#2" 

PRINT*," 3-DON' T STORK" 

READ*,QUEB 

USE SUBROUTINE TRNSPLT TO SWITCH TOE AVERAGE PLOT WTO A PLOT BAY. 

IF (QUES .EQ. 1) CALL IRN8PLT(MACAVB, FREQ, AVPARMS , SAVMACl , SAVFRQl, 

* PAKAHS1) 

IF (QUES .EQ. 2) CALL TRNSPLT(HACAVE, FREQ, AVPARMS ,SAVMAG2 .SAVFRQ2 , 

* PAHA MS 2) 

LOO CALL PRTS TAT (PAR ANSI , PARA MS 2 , WKPARHS ) 

PRINT*, '1)0 YOU WISH TO CONTINUE' TO PROCESSING? 1-YES " 

PRINT*," 2-NO, GO BACK 

PRINT*," POR MORE SPECTRA 

PRINT*," 3-QUIT SPECTRA" 

READ*, QUES 

IF (QUES .EQ. 2 )COTO HO 
17 (QUES .EQ. 3) STOP 

BEGIN PROCESSING 

ASK THE USER WHICU PROCESS IS DESIRED. THE SUBROUTINE MENU 
ASKS THIS AND RETURNS WITH A NUMBER IN QUES. 

L05 CALL MENU (QUES) 

ANY SUBROUTINE WITH "BET" IN IT MERELY SETS UP A CALI. TO TOE 
ACTAUL MATHEMATICAL ROUTINE DESIRED. 

IFfQUEB .EO. 1)CALL DIFF8ET(SAVHAG1 .SAVFRQl ,PARAMB1 .BAVMAG2 , 

* SAVFRQ2 , PARAMS 2 , KOftKMAG , WORKFRQ , WKPARHS ) 

IP (QUES .EQ. 2 )CALL WND08ET(SAVMAG1., SAVFRQl , PARA M3 1 .8AVHAG2 , 

* 8AVFRQ2 .PARAMS2 , 

* WORKMAG, WORKFRQ, WKPARHS) 

IF(QUES .EQ. 3) CALL TENS SET (SAVMACl .SAVFRQl .FARAHSl .8AVMAC2,, 

* SAVFRQ2 ,?ARAHS2 , 

* WORKMAC , WORKFRQ , WKPARHS ) 

IF(QUE8 .EQ. A)CALL AVRGSF.T ( SAVMACl , SAVFRQl .PARAHS1 ,SLTJJJS2 , 

* SAVFRQ2.PARAMS2, WORKMAG, WORKFRQ, WKPARHS) 


IP(QIIES .EQ. 5)CALL PRTS TAT ( PARAMS 1 .PARAMS 2, WKPARMS) 


FFT8 IZE-IPIX(WKPARMS (3 ) ) 
PLTSIZE-IFIX(WKPARM8(2) ) 
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IF(QUEB .EC). 6) CALL PLOTTER { WORKMAN , WORKFRQ , FFTS I ZE , 

* PLTSIZE.HAXHAG.HINMAC.HAXFREQ, 

* OPTION .SLOPE ,YINT,FIT1 ,LOLIMIT,HILIMIT, 

* WKPARMS (1 ) , 8H0TNUH, ALTITUD .ROCKOPD , IUR6 , 

* IMIN ,6EC8 , DIGRATE f OVLP , GA’N , JJ ) 


IP(QUE8 .EQ. 8) CALL C0PY8ET(SAVMACl .BAVFRQl .PARAMS1 , 

* SAVMAG2 ,SAVFRQ2 .PARAMB2 , 

* VORKMAG, WORKFRQ .WKPARMS) 

IF(QUE8 .CT. 8)CALL PRNTERRO) 

C 

C IF QUES WAS NOT EQUAL TO 7, CO BACK TO THE MENU. 

C 

IF (QUE8 .SE. 7) GOTO 105 


C 

C GO BACK TO THE BEGINNING OF SPECTRA FOR MORE PLOTS 
C 

110 PRINT*, "ARE YOU FINISHED? 1-YES" 

PRINT*," 2-HO" 

READ*, QUES 

IF (QUES .EQ. 1) STOP 

PRINT*, "GO BACK TO MAKE MORE PLOTS." 

PRINT*, ,i DIFFERBJT TIME, SAME ANALYSIS? 1-YE6" 

PRINT*," 2-NO" 

READ*, QUES 
REPEAT* 0 

IF (QUES .EQ. 1) REPEAT-1 
GOTO 10 


STOP 

END 

C 

C SUBROUTINE LOAD READS A SPECIFIED CHANNEL ON A 


C 

c 

C 

c 

c 

C 

c 

C 

C 

c 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PDP2CSC CONVERTED TAPE AND STORES THIS DATA IN ARRAY 
DATA. THE NUMBER OF DATA POINTS TO BE READ AND THE 
TIME CORRESPONDING TO THE FIRST DATA POINT ARB ENTERED 
BY THE USER. IF THE STARTING TIME ENTERED DOES HOT 
CORRESPOND TO A SPECIFIC DATA POINT, THE FIRST DATA 
POINT OF THE DESIRED CHANNEL FOLLOWING THE GIVEN 
STARTING TIME WILL BE THE FIRST TO BE READ. FOLLOWING 
IS A LIST OF THE SUBROUTINE VARIABLES. . 

(NOTE — THE NUMBER IN PAR BJ THESES IN THE READ STATEMENT 
SPECIFIES THE DEVICE TO BE READ. THIS SHOULD BE WDIFIED 
IN ACCORDANCE WITH THE UBER'B PROGRAM STATEMENT.) 

DATA -ARRAY CONTAINING THE CHANNEL SAMPLES. 

N .-NUMBER OF SAMPLES TO BE READ. 

LTIHE -LAUNCH TIME IN UNIVERSAL TIME (IN SECONDS). 

STIME -TIME CORRESPONDING TO THE FIRST SAMPLE TO BE 
READ (IN SECONDS FROM LAUNCH). 

PTIHE -TIME CORRESPONDING TO FIRST SAMPLE IN A RECORD. 

8MPTIHE -CHANNEL SAHPLINC PERIOD. 

TIME -TIME OP DATA POINT IN UNIVERSAL TIME (LTIHB+ 

STIME ) . 

FRAME -THE FRAME NUMBER (A FRAME IS A SET OF FIVE DATA POINTS). 

CHANNEL -THE DESIRED CHANNEL TO BE READ. 

POINT -THE DATA POINT CORRESPONDING TO TIME 

IPOINT -THE LARGEST INTEGER LESS THAN OR EQUAL TO POINT. 

STPOINT -THE FIRST DATA POINT TO BE READ (STARTING POINT). 


SUBROUTINE LOAD (DATA.N, CHANNEL, 1TIME, STIME, SMPTIME) 

REAL DATA (N ) .ARRAY (1 251 ) , LTIME , STIME ,T ME , 8MPI IME , PT IME .POINT , 
* TIKE LI 

INTECER IARRAYU250) .IPOINT, STPOINT, CHANNEL 
EQUIVALEN CE ( ARRAY ( 2 ) , IARRAY ( 1 ) ) 

TIME-LTIME+STIME 
READ(3) ARRAY 


PTIME-ARRAY(l) 

READ<3) ARRAY 

BMP! IKE“(ARRAY(1)-PTIHE ) /250 ,000000 
REWIND 3 
C 

C LOCATE THE PROPER RECORD 
C 

10 READ(3) ARRAY 

C T1HEL1 IS THE TIME OF THE RECORD PREVIOUS TO THE DESIRED ONE. 

C THE .0001 IS USED TO PREVQ1T ROUND OFF ERRORS. 

TIMEL1 - TIHE-(SHPTIME*250.0000) ♦ .0001 
IF CARRAY(l) .CT. TIMEL1) GOTO 20 
GOTO 10 

20 POINT- (TIHE-ARRAY(l))/(BMPriME/5.0)*1.0 
IPOINT-POINT 

IF(POINT-FLOAT(IPOINT) .LT. 0.000001) GOTO 30 

IPOINT-IPOINT+1 

IFdPOIKT .LE. 1250) GOTO 30 

READ<3) ARRAY 

IPOINT-l 

30 FRAME- (IPOINT-l )/5+l 

STPOINT»(FRAHE-l)*5+CHANNEL 
IF(IP01NT .LE. STPOINT) GOTO 40 
STPOINT-STPOINT+5 
1F(ST?0INT .LE. 1250) GOTO 40 
READ (3) ARRAY 
STPOINT-CHANNEL 
C 

C READ IN THE DATA 
C 

40 DO 100 I-l.N 

DATA(l)-FLOAT (IAR3AY ( STPOINT ) ) 

8TPOINT-STPO INT+ 5 
IF{ STPOINT .LE. 1250) GOTO IOC 
READ<3) ARRAY 
STPOINT-CHANNEL 
100 CONTINUE 
RETURN 
END 
C 
C 
C 
C 
C 

C SUBROUTINE CALIB PERFORMS THE NECESSARY 
C CALIBRATIONS SO THAT THE DATA IS REPRESENTATIVE 
C OF THE OUTPUT FROM THE EXPERIMQ1T. THESE 
C CALIBRATIONS INCLUDE ADJUSTMENTS FOR EXPERIMENT 
C GAIN, DIGITIZER GAIN, SIGNAL OFFSETS, AND 
C DISCRIMINATOR NON-LINEARITIES. THE NON-LINEARITIES 
C ARE REMOVED WITH A CALIBRATION CURVE OBTAINED BY 
C A PIECEW IB E-LIN EAR FIT TO THE CALIBRATION VALUES 
C (BANDCAL(5) ). 

C 

C 

SUBROUTINE CALIB (DATA, N ,BANDCAL, GAIN) 

REAL DATA(N) , BAN DCA L( 5 ) , GAIN 
INTEGER N 
DO 100 1-1 ,N 

IF(DATA(I) .GT. BANDCAL(2)) GOTO 10 
DATA(l)-((l .25/(BANDCAL(2)-BANDGAL(l) ))*DATA(I) 

* -(1.25*BANDCAL(l)/(BANDCAL(2)-BANDCAL(l)))-2.5)/GAIN 
GOTO 100 

10 IF(DATA(I) .CT. BANDCAL<3>) GOTO 20 

nATA<I)-({ l . 25/ (BANDCALC3 )-BANDCAL(2'j ) )*DATA(I ) 

* -(1 .25*3ANDCAL(2)/(BANDCAL(3)-BANDCAL(2) ))’-l .25) /GAIN 
GOTO 100 

20 IF'.DATA(I) .GT. BANDCAL(4)) GOTO 30 ’ 

DATA(l)-((1.25/(BANDCAL(4)-BANDCAL(3)))*DATAd) 

* -(1.25*BANDCAL(3)/(BANDCAL(4)-BANDCAL(3))))/GAIN 
GOTO 100 

30 DATA(I)-((1 .25/ (BANDCAL(5)-BANDCAL(4)))*DATA(I) 

* -(1 .25*BANDCAL(4)/(BANDCAL(5)-BANDCAL(4)))+l .25) /GAIN 


100 CONTINUE 
RETURN 
END 
C 
C 

c 

c 

c 

C SUBROUTINE SpECTRM SEGHENTS THE DATA, REMOVES TOE MEAN, 

C WINDOWS TOE DATA, COMPUTES TOE PERIODOGRAM, UPDATES AN 
C AVERAGE OP THE PERIODOGRAH8 , AND THEN DISPLAYS TOE 

C INDIVIDUAL PERIODOGRAHS . AN AVERAGED PERIODOGRAM IS 

C OBTAINED UPON RETURN FROM THE SUBROUTINE. 

C 

C 

SUBROUTINE S PECTRM (DATA , MAG, MAGAVE , PREfl ,X,Z,N,K,H,M2, L ,OVLP , 

* WINDONO , SMPTIME .OPTION .DISPLAY .NUMBER, HINHAG, 

* PIT2 , PIT1 , LOLIHIT , HILIMIT , FC , FAC , FLE , BEAN D , 

* LT IME , ST IHE , SHOTO UM, DIGRATE, TOTAL .PLOTS ) 
INTEGER N ,K,M,H2 ,OVLP,D, L, WINDONO ,IWK(1 2) .DISPLAY .OPTION , 

* NUMBER , S .BRAND , V , FIT1 , FIT2 .PLOTS ( 100 ) .TOTAL , PRES PLT , DIGRATE 
REAL DATA(N) ,HAG(M2) ,MAGAVE(M2) ,FREQ(M2) ,X(M) .MEAN ,WK(1) , 

* SHPT IME .MAXMAG, MINMAG , MAXFREQ , LOLIMIT , HILIMIT , FC , FAC , FLE , 

* LTIKS.STIME, SLOPE, YINT 
COHPLEX Z(M2) 

C 

C INITIALIZE ARRAYS 
C 

DO 100 I-1.M2 
MAGAVE (I)“0 .0 
Z(l)-((0.0,0.0)) 

100 CONTINUE 

DO 700 I-I.H 
X(l)«0 .0 
700 CONTINUE 

D-(N-OVLP)/K 

S-l 

PRES PLT- 1 
C 
C 

C BEGIN SPECTRAL ESTIMATION 
C 

PRINT* .PLOTS Cl), PLOTS ( 2) .PLOTS (3 ) , PLOTS (100) 

C 

DO 200 J-l.K 
C 

C THE ARRAY PLOTS CONTAINS THE PLOT NUMBERS DESIRED BY THE USER. 

C THE FOLLOWING STEP CHECKS WHETHER TOE DxSIRED PLOT NUMBER (COH- 
C TAINED IN PLOTS (PRES PLT) ) MATCHES TOE J'TO PLOT ABOUT TO 
C BE CALCULATED BY THE "DO 200" LOOP. IF THERE IS A HATCH THE 

C LOOP IS EXECUTED AND PREBPLT IS INCREMENTED. IF THERE IS NO 

C MATCH THE LOOP IS SKIPPED ( J BEING AUTOMATICALLY INCREMENTED) 

C AND AGAIN THE STEP SEES IF PLOTS (PRES PLT) ■ J. 

C A MODIFICATION BY ROD STOLTZFUS, 2/16/84 
C 

FRINT*. PLOTS (PRES PLT) , J 

IF (PLOli(PRESPLT) .NE. J) GOTO 200 

PRESPLT-PRESPLT+1 

C 

C READ IN TOE APPROPRIATE DATA SEGMENT. 

C 

DO 300 I-l.L 

X(I)“DATA(I+(J-1)*D) 

300 CONTINUE 

C 

C MEAN REMOVAL 
C 

MEAN-0.0 
DO 400 I-l.L 

HEAN-MEAN+(X(I)-HEAN) /FLOAT(I) 
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400 CONTINUE 

DO 500 I-1,L 

x(i)-x(i)-hkah 

500 CONTINUE 
C 

C APPLY WINDOW TO DATA. 

C 

CALL WINDOK(X,L.M,VXNDOno) 

C 

C CONFUTE PERIODOCRAM 
C 

C FOR MORE INFORMATION ON PITRC 8EE THE IMBL DOCUMENTATION. 

C 

CALL FFTRC(X P M,Z,IWK,WK) 

C DO 800 I- 1, M2 (OPTIONAL) 

C Z(I)-C0NJC(Z(I)) (OPTIONAL) 

C800 CONTINUE (OPTIONAL) 

DO 600 1-1 ,M2 

HAC(I)-(2.0*8MPriKE/FLOAT(L)J*(CABB(Z(l)))**2 
FREQ( I)-(I-l) /(FLOAT (H)*SHPTIME ) 

600 CONTINUE 
C 

C AVERACE PERIODOGRAMS 

C ] 

IK) 900 I-1.M2 

MAGAVE ( 1 ) ■ (1 . 0/ FLOAT (TOTAL ) )*HAC( I )+MACAVE (I) 

900 CONTINUE 
C 

C DISPLAY INDIVIDUAL PERIODOGRAM 
C 

IF(DISPLAY .EQ. 2)GOTO 200 

IF(S .GT. NUMBER) GOTO 200 

CALL FILTCAL(HAG,FREQ,M2,FC,FAC,PLB) 

CALL DEEMPH(MAC,FREQ,M2,BBAND) 

CALL HAXMIN(HAG,FHEQ,MAXMAG,MAXFREQ,H,H2) 

CALL PLOTTEROIAG, FREQ.M, M2 ,HAXHAC,HINHAG,HAXFREQ, OPTION, SLOPE, 

* YINT.FIT2 ,LOLIMIT,HILIMIT,8TIKE .BHOTNUH.ALTITUD, 

* ROCKS PD, IHRS, IMIN , SECS , DIGRATE ,OVLP , GAIN , J) 

IF(FIT2 .EQ. 2)GOTO 10 

PRINT*, "DO YOU WISH TO STORE THE SPECTRAL INDEX7 1— YES" 

PRINT*," 2— NO" 

READ*,V 

IF( V .EQ. 1)CALL STORE (SLOPE , YINT, LTIME ,STIME , SHPTIME , N) 

10 S-S+l 
200 CONTINUE 
C 

C RECONSTRUCT ARRAY FREQ( ) SINCE SUBROUTINE PLOTTER DESTROYS IT. 

C FOLLOWING LOOP ADDED ON B. TOMEI'S SUGGESTION. 10/24/83 

C 

DO 999 I-1.M2 

FREQ( I )«(I-1) / (FLOAT(M)*SMFTIHE ) 

999 CONTINUE 1 

RETURN 3 

END 

C | 

C | 

c 

C 

c 

C SUBROUTINE KAXHIN CALCULATES THE VALDES OF HA3SHAG, MAXMIN, 

C AND MAXFREQ AND ROUNDS THEM TO APPROPRIATE VALUES FOR AXES 
C PLOTTING. 

C ' 1 

C | 

SUBROUTINE MAXMIN (HAG, FREQ ,MAXHAG, MAXFREQ ,M ,H2 ) ij 

INTEGER M,M2 1 

REAL MAG(K2),FREQ(M2),MAXMAG, MAXFREQ i 

MAXMAG-0 .0 I 

MAXFREQ-FREQOI2) . 3 

DO 100 I-1.M2 | 

IF(HAG(I) ,GT. MAXMAG) MAXMAG-MAG( I ) 3 
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'MW' " 
A- '.Ss., 


c 

c 

c 

c 

c 

c 

c 

c 

c 

0 

c 

c 

c 

0 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


100 CONTINUE 

CALL LIMITS (HAXHAC) 

RETURN 

END 


, SUBROUTINE FILTCAL PERFORMS A CALIBRATION TO THE SPECTRUM 
IN ORDER TO COMPENSATE FOR THE ROLLOFF EFFECTS OF THE 
TELEMETRY FILTER, TOE AC AMPLIFIER, AND THE LOGARITHMIC 
ELECTROMETER. THE TELEMETRY FILTER IS MODELED AS AN 
IDEAL 3-POLE LOW PASS FILTER (SEE ZIMMERMAN ,1980) , THE 
AC AMPLIFIER AND LOG ELECTROMETER ROLLOFFS ARE MODELLED 
AS SINGLE POLE LOW PASS FILTERS. THE INPUT PARAMETERS 
ARE 

FC — TOE UPPER 3DB FREQUENCY OF THE TELEMETRY FILTER. 

FAC — TOE UPPER 3DB FREQUENCY OF THE AC AMPLIFIER. 

PEL — THE UPPER 3DB FREQUENCY OF THE LOG ELECTROMETER. 

HC,HAC,HEL ARE THE MAGNITUDES SQUARED OF THE TRANSFER FUNCTIONS 

CORRESPONDING TO WE RESPECTIVE FILTER. 


SUBROUTINE FILTCAL(HAG,FREQ,H2 ,FC,FAC ( FLE) 
INTEGER M2 

REAL MAG(H2), FREQ( M2), OMEGA, PJ.WC.WAC.WLE 
Fl-3. 141 59265 
WC«2.0*PI*FC 
WAC-2 ,0*PI*FAC 
WLE-2.0*PI*FLE 
DO 100 I-1,M2 
OHECA-2.0*P1*FREQ(I) 

HC- 7 . 7 53 38 *WC**6 / ( ( ( 2 . 10 19 9*WC**2-OHECA**2 )**2 

* +4,40165*<WC*OHECA)**2>*<1 ,754B0*VC**2 

• +OMEGA**2) ) 

HAC-WAC**2/ <WAC**2+OKEGA**Z) 

HLE»WLE**2/ (V!LE**2+0MEGA**2) 

MAG ( I ) -KAG( I ) / (HCTOAC *HLE ) 

100 CONTINUE 
RETURN 
END 


SUBROUTINE DEEMPH PERFORMS DE- EMPHASIS ON THE SPECTRUM 
OBTAINED PROM TOE FINE STRUCTURE EXPERIMENT. TOE PRE-EMPHASIS 
FILTER 1$ MODELED AS A SINGLE-POLE HIGH PASS FILTER WITH 
CUTOFF FREQUQICY EQUAL TO THAT OF THE BROADBAND AMPLIFIER. 

TOE DE-EMPHASIS FILTER 18 THEN OBTAINED BY TAKING THE INVERSE 
OF THE PRE-EMPHASIS FILTER TRANSFER FUNCTION. DE-EMPHASIS 
IS REQUIRED FOR PAYLOADS 31.028 , 31.029, AND ANY OTHER 
PAYLOADS WITH TOE SAME EXPERIMENT CONFIGURATION, THE 
PARAMETERS ARE 

1) BRAND — THIS NUMBER REFERS TO TOE NUMBER OF TOE 

BROADBAND AMPLIFIER WHOSE OUTPUT IS DNDER 
ANALYSIS. FOR BRAND .Eq. 0, DE-Q1PHASIS 
IS BYPASSED. * 

2) WC — THIS IS THE LOWER 3DB CUTOFF FREQUENCY OF THE 

BROADBAND AMPLIFIER. 

3) HDE — THIS IS THE MAGNITUDE SQUARED OF THE DE-EMPHASIS 

FILTER TRANSFER FUNCTION. 


SUBROUTINE DEEMPH (MAG,FREQ, M2 ,BBAND.) 
INTEGER M2 , BRAND, L 

REAL MAG(M2),FREQ(M2), PI, OMEGA, WC, HDE 

PI-3.14159265 

IP(BBAND .EQ. OJREIURN 
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c 

C PERFORM DE- EMPHASIS 

c 

WC-2.0*P1*400.0 

IF( BRAND .EQ. 2) WC"2.0*PI*53 .9 
DO 100 1-2, M2 
OMBCA-2.0*PI*PREQ(l) 

HDE-(OMECA**2+VC**2) /CMEGA**2 
HAG(I)-HDE*MAG(l) 

100 CONTINUE 
RETURN 
END 
c 
c 
c 
c 
c 

C SUBROUTINE PLOTTER PLOTS THE POWER DENSITY VERSUS 

C FREQUENCY OR SCALE SIZE IN LINEAR OR LOGAR1THHIC FORM. 

C A LINEAR FIT MAY ALSO BE PERFORMED ON THE LOGARITHMIC PLOT. 

C 

(£**** 

C NOTE THAT THIS SUBROUTINE WILL ROUND OFF POINTS WHICH ARE OFF 
C THE PLOTTING SURFACE. THUS THE ROUTINE IS A BIT DESTRUCTIVE. 

C A LOCAL COPY OF THE MAGNITUDE ARRAY COULD BF, BE MADE IF ABSOLUTE 
C INTEGRITY OF THE DATA IS REQUIRED. 

C 

C IF CURVE FITTING IS CHOSEN, THE INPUT MAGNITUDE INFORMATION IS 
C TOTALLY DESTROYED IN PROVIDING INFORMATION TO THE LINE FITTING 
C SUBROUTINE. 

C***« 

C 

c 

c 

SUBROUTINE PLOTTER (HAG,FREQ,M,M2,MAXHAG,MINMAC,HAXFREQ ( OPTION, 

* SLOPE , YINT , FIT ,LOL IMIT , HILIMIT , BT IME , SHOTOUH , 

* ALTITUD, ROCKS PD, IHRS ,IMIN,SEC8 ,DIGRATF,,OVLP, 

* GAIN ,FLOTNUM) 

INTEGER M ,H2 , OPT ION , FIT , J , IHRS , IMIN , 0 VLP , DIGRATE , PLOINUM 
REAL HAG( M2 ) , FREQ(M2 ) , MAXMAG , MINMAG , MAXFREQ , YINT , SLOPE , 

* LOLIHIT ,HILIMIT , Xl < 3 ) ,Y1(3 ) ,STIME .TIME .ALTITUD , aOCKSPD , 

* SECS , GAIN . SHOTN UM.REALK , REALM , RELOVLP , REALDIG , LOCMAG( 20 A 9 ) 

C 

C 

C 

c**** 

C SUBROUTINE PLOTTER ADJUSTS THE DATA TO FIT INTO THE LEGAL BOUNDS 
C OF THE REQUESTED PLOT. FOR EXAMPLE, IT CHANGES ANYTHINC < 

C MINMAG TO MINMAG. TO PRESERVE THE ORIGINAL DATA THE ORIGIN/ L 
C DATA IS SAVED IN HAG. 

C 

DO 5 I-1.M2 
LOCMAG(l)-MAGd) 

5 CONTINUE 
C 

C INITIALIZE PLOTTING PARAMETERS 

C 

C 

CALL U8TART 

C THE FOLLOWING DISPLAY CAPABILITY WAS ADDED BY RODNEY 

C STOLTZPU8. 10/09/83. 

C TO USE, REMOVE THE C'S AT THE BEGINNING OF THE LINES. 

C ALSO MAKE SURE APPROPRIATE PARAMETERS ARE PA8BF.D INTO SUBROUTINE 

C PLOTTER AND SPCTRH. 

CALL UWRITE(52 . ,98 . , "ROCKET BHOT i "> 

C CALL UVRITE ( 52 . , 9 5 . , "LAUNCH TIME ;") 

CALL UWRITE C 52 . , 9 5 . , "PLOT NUMBER;") 

CALL UWRITE (52., 92., "SEC AFTER LAUNCH;") 

C CALL UWRITE(52., 89., "SPEED ;") 

C CALL UWRITE(52., 86., "ALTITUDE ;") 

C CALL UWRITEC52., 83., "ROCKET ANGLE ;") 

CALL UWRITE (52 . ,80 . , "FFT BLOCK LENGTH;") 

CALL UWRITE ( 52 . , 7 7 . , "SEGMENT OVERLAP;") 
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CALL UWRITE(52., 75. ."DIGITIZING RATE ; ") 

CALL UWRITEC52., 72., "ROCKET SPIN RATE;") 

CALL «HRITE{52., 69., "GAIK ;") 

aLL UPBET( "PRECISION" ,6 .) 

CALL USETC'REAL") 

CALL UWRITE (8 4 . , 98 . , SHOW UM ) 

CALL UWRITE(84.,72.,SPINRAT) 

CALL UWKITE(84 , ,B3 , , ANGLE) 

CALL UWRITE(84.,B6.,ALTITUD) 

CALL UWRITE (64., 89. , ROCKS PD) 

STIHE IS THE TIME AFTER LAUNCH 
WHEN STIHE IS PASSED TO SUBROUTINE PLOTTER IT CONTAINS THE 
TIME OF THE FIRST PLOT. IN PLOTTER, A LOCAL COPY CALLED TIHE 
IS USED FOR LOCAL HANDLING, 

IF PLOTNUH IS 0 THE AVERAGE PLOT IS BEING REQUESTED. THUS 
CORRECT TIHE S.T. IT IS THE FIRST PLOT'S TIHE (I.E. STIHE). 

THE FOLLOWING FORHULA ASSUMES FFT LENGTH “ SEGMENT LENGTH, 

TIHE“STIKE+(( FLOAT (PLOTNUH- 1 ) ) * < FLOAT (M-OVLP ) ) / (FLOAT (DI GRATE ) ) ) 
IF( PLOTNUH .EQ, 0) TIME-BTIHE 


CALL UWRITE(84.,92.,TIHE) 
aLL UWRITE (96 . , 95 . , SECS ) 
aLL UWRITE(84.,69. , GAIN) 
aLL USBT(" INTEGER") 

MUST MAKE ALL VARIABLES REAL TO BE PASSED WHEN PLOTTING. 


t 

( 

t 

} 


aLL UWRITE(84 . , 95 . , FLOAT (PLOTNUH) ) 
aLL UWRITE(84 . , 95 . , FLOAT ( IHRS ) ) 
aLL UWRITE (90., 95., FLOAT (IMIN)) 
aLL UWRI;:e ( 84 . , 72 . , FLOAT (DIGRATE ) ) 
aLL UWRITE(84. ,77, , FLOAT (OVLP) ) 
aLL UWRITE (84., 80., FLOAT (H)) 

CALL UBET("iEXT") 

END DISPLAY MODIFiaTIONS BY R.S. 


aLL USET("XBOTHLABEL") 
aLL U8ET("YB0THLABEL") 
aLL USETC'OWNSaLE") 
aLL UPSET("XLABEL"," 5 ") 
aLL UPSET ("YLABEL"," ;") 

CALL UPSETCTERHINATOR”,";") 
aLL UPSET("VERTiaL",3.0) 

CALL UPSET("HORIZONTAL",2 .0) 

ai.L ubet("sofiware") 

aLL UDAREA(0.0,7.0,0.0,5.7) 

aLL UPRINT(40 .0 ,2 .0 , "FREQUENCY (HZ) ; ") 

aLL UM3VE(4. 0,0.0) 

aLL UROTAT(90 .0) 

aLL UPRINT( 25. 0,0.0, "POWER DENBITY (DELTA N/N) /HZ;") 
aLL UFRINT(79.0,1.5,"2;") 
aLL USET( "HARDWARE") 

CALL UBETC’SYSTEH") 

THE FOLLOWING THREE STATEMENTS MAY BE USED TO 
SUPERCEDE THE OUTPUT OF SUBROUTINE HAXMI1I. IF 
THIS IS DESIRED, REMOVE THE COMMENT MARKER8. 

MINHAG-10E-12 (OPTIOHA7.) 

HAXMAG“10E-4 

HAXFREQ-1000.0 (OPTIONAL) 

aLL UWINDO(0.0,MAXFREQ,0.0 ,MAXHaG) 

IP(0PTI0H .EQ. 2)GOTO 10 

PERFORM LINEAR PLOTTING 


CALL UFL0T1 ( FREQ , MAC, PLOAT<M2 ) ) 

GOTO 20 
C 

C PERFORM LOGARITHMIC PLOTTING 
C 

10 CALL UWINDoO .0 .MAXFREQ .MINHAG.MAXHAG) 

CALL USET ( "XYLOGAR ITHMIC " ) 

CALL US ET ( "LOGXAXIS " ) 

CALL USETO'LOGYAXIS") 

J-0 

DO 100 1-1 ,H2 

IP(MAC(I) .LE. MIN HAG) MAG( I ) -MINHAG 
IF(FREQ(I) .LE. MAXFREQ) J-J+l 
100 CONTINUE 

FREQCD-0.01 

CALL UPLOTl(FREQ,MAG,FLOAT(J)) 

CALL UFLUSH 
CALL UPAUSE 
C 

C RESTORE THE HAG USING THE SAVED LOCHAC, 

C 

DO 195 1-1, HZ 
MAC(I)-LOCMAGd) 

195 CONTINUE 
C 
C 

C PERFORM LINEAR LEAST-SQUARES FIT IF FIT .NE. 2. 

C 

c 

IP(FIT .EQ. 2) RETURN 
FREQ(1)-0.0 
C 

C FLAG DATA POINTS WITH MAGNITUDE EQUAL TO EERO 
C OR WITH FREQUENCY OUT OF LIMITS . 

C 

C **** ORIGINALLY B. TOMEI HAD .EQ. IN THE SECOND PAIR OF FOLLOWING 
C ****** CONDITIONALS. THEY HAVE BEEN CHANGED TO .LE. 

C***** R.s. 3/11/84 
C 

DO ZOO 1-2 ,M2 

IF(FREQd) .LT. LOL IMIT ) FREQ ( I ) -0 . 0 
IF(FREQd) .CT. HIL1MIT )FREQ( I ) “0 .0 
IF(HAGd) .LE. HINHAG)FREQd)-0.0 
IF(MAG(I) .LE. 0.0)FREQ(I)-0.0 
200 CONTINUE 
C 

C REMOVE FLAGGED DATA POINTS AND COMPRESS DATA 
C ****FLACGED POINTS ARE THOSE WITH FREQUENCY LESS 
C ****THAN 0.1. IF THE PLOT'S LOWEST FREQUENCY IS 

C ****LESS THAN 0 .1 CHANGE THE CONDITIONAL WHICH 

C **** WHICH REMOVES THE DATA POINT 

C 
C 

.1-1 

DO 300 I-1.H2 

IF(FREQ(I) .LE. 0.0)C0TO 300 
MAG( J)-ALOGlO(MACd)) 

FREQ(J)-ALOG10(FREQ(I)) 

J-J+l 

300 CONTINUE 
C 

C PERFORM LINEAR FIT 
C 

CALL ULIKFTCFREQ.MAG.FLOATCJ-l) .SLOPE, YIN1 

c 

C GENERATE FITTED CURVE (THREE POINTS ONLY) 

C 

XI ( 1 ) "LOLIHIT 

Y1 (I ) - ( 1 0 ,0**YINT )* (XI (1 ) **SLOPE ) 

XI ( 2 ) - ( LOLIMIT+BILIHIT ) / 2 .0 

Yl ( 2 > -(10 ,0**YINT )*( XI ( 2 )*«SLOPE ) 

X1(3)-HILIMIT 

Y1(3)-(10.0**YINT)*(X1(3)**SLOPE) 
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C PLOT FITTED LINE 
C 

CALL USET("LINE") 

CALL USETC'OLDBCALE") 

CALL GPL0T1 (X).,Y1 ,3 .0) 

20 CALL UFLUSH 
CALL UPAUSE 
CALL UERASE 

RESTORE ARRAY MAC USING THE SAVED DATA IN LOCKAC 

DO 25 I-1,H2 
MAC(l)»LOCHAC(I) 

25 CONTINUE 

CALL VEND 

RETURN 

END 


SUBROUTINE STORE BT0RE8 THE SPECTRAL INDEX DATA IN ARRAY 
INDEX IN THE FOLLOWING FORMAT- TIME .SPECTRAL INDEX, 
MAGNITUDE. ARRAY INDEX THEN SERVES AS THE INPUT ARRAY 
FOR PROGRAM SPINDEX. 


SUBROUTINE STORE (SLOPE, YINT, LTIME .BTIHE.SHPTIME.N) 
INTEGER N 

REAL SLOPE , YINT, LTIME , STIME , SHPTIME , AVTIME 
AVTIME- (LTIHE+STIHE )+(FLOAT(N)*BMPT IME ) /2 .0 
WRITE(4,11) AVTIME .SLOPE , ( 10 .0**YINT ) 

II F0RMAT(1X,F9.2 , IX, FlO. 5 , IX, El 3. 6) 

RETURN 

END 


SUBROUTINE WINDOW APPLIES A SPECIFIED WINDOW TO THE DATA 
SEGMENT X, THE WINDOW 18 SPECIFIED BY WINDONO. THE 
WINDOWED DATA IS ALSO NORMALIZED BY THE INCOHERQtT 
POWER GAIN OF THE WINDOW (IPG). 


SUBROUTINE WINDOW (X,L,M, WINDONO) 
INTEGER L,J,R, WINDONO, H 
REAL X(M) ,PI ,IPC 
PI-3.141592653 

GOT0(10 ,20,30 ,40 ,50 ,60)WINDON0 

RECTANGULAR WINDOW 

10 RETURN 

BARTLETT WINDOW 

20 J-L/2 
R-J+l 
IPG-0.0 
DO 100 1-1 , J 

WEIGHT-FLOAT ( 1-1 ) /FLOAT (J) 
X(I)»X(I)*WEIGHT 
IPG-WEIGHT**2+IPG 
100 CONTINUE 

DO 200 I-R,L ■ 

WEIGHT- FLOAT(L- I+l ) /PLOAT(J) 

X(I)-X(I)*WEIGHT 

IPG-WEIGHT**2+IPG 
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200 CONTINUE 
IPC-IPC/fLOAT<L) 

DO 201 1-1,1 

X(I)-K(I)/8QKT{IPC) 

201 CONTINUE 
RETURN 

C 

c hanning window 

c 

30 IPG-0.375 
DO 300 I-l.L 

WE1CHT-0 . 5* ( 1 .O-COS(2.0*PI*(l-l)/L) ) 
X(I)-X(I)*WEIGHT/BQRT(IPC) 

300 CONTINUE 
RETURN 
C 

C HAMMING WINDOW 
C 

40 IPG-0.3974 

DO 400 I-l.L 

WEICHT-0 .54-0.46*C0S(2,0*PI*(I-l) /L) 
X(I)-X(I)*WEICHT/SQRT(IPG> 

400 CONTINUE 
RETURN 
C 

C BLACKMAN WINDOW 
C 

50 IPG-0.3046 

DO 500 1-1 ,L 

WEIGHT-0. 42-0. 5*COB(2.0*PI*(l-l)/L)+0.08*COS(4.0*Pl*(I-l)/L) 

x(i)-x(i)*weight/sqrt(ipg) 

500 CONTINUE 
RETURN 
C 

C BLACKHAN-HAKRIS WINDOW (4-TERM) 

c 

60 IPG-0 .269672277 
DO 600 I-l.L 

WEIGHT-0 .40217-0.4970 3*C0S ( 2 .0*P I*FLGAT( 1-1 ) ft) 

* +0 .09392*008(4 .0*PI*FL0AT (I-l)/L) 

* -0 .001B3*C08(6 .0*PI*PLOAT(I-1) ID 
X(I)-X(I)*WEIGHT/8QRT(IPG) 

600 CONTINUE 
RETURN 
END 
C 
G 
C 
C 

c 

C SUBROUTINE LIMITS ROUNDS THE VALUE OP MAX FOR EASE IN AXES 

C PLOTTING. FOR 10**(N+1) > MAX > 10**N ,(N GREATER THAN OR 

C EQUAL TO 1.0), LIMITS ROUNDS THE VALUE OF MAX U? TO THE 
C NEXT (N-l)TH POWER OF TEN. FOR 10**-(N-1) > MAX > 10**-U , 

C (II LESS THAN 1.0), LIMITS ROUNDS THE VALUE OF MAX UP TO 

C THE HEX! -(N+1)TH POWER OF TEN. 

C 

C 

SUBROUTINE LIMITS (MAX) 

REAL Z,Y, FACTOR, MAX, DELTA 
IF (MAX ,GE. 10.) GOTO 10 
DO 100 1-1,15 
Z-10 .0**(I-14) 

100 IF (MAX .GT. Z)FACT0R-10 .0**(I-14) 

Y- (AINT( 10 . 0*MAX/FACT0R ) ) I 10 . 0 
MAX-FACTOR* (Y+.l) 

RETURN 

10 DO 200 1-1,6 
Z-10 .0**1 

200 IF (MAX .GT. Z)FACTOR»10 .0**(I-1) 

DELTA-0.0 

IF(AMOD(MAX, FACTOR) .GE. 0 .000.1 )DELTA-1 .0 


on o o o o non o o n o o o oo o o o 


>IAX-FACTOR*(AINT(MAX/FACTOK)*DELTA) 

RETURN 

END 

C 

C 

c 

c 

c 

c 

c SUBROUTINE SELECT ASKS FOR THE PLOT NUMBERS DESIRED. IF THE 

C USER ONLY HANTS TO SEE PLOTS 3 AND 12, TOE USER BITERS THESE 

C NUMBERS. TOE PLOT NUMBERS MUBT BE BiTERED IN ORDER, 1,2. PLOT 

C 12 CANNOT BE REQUESTED BEFORE PLOT 3. TOE SUBROUTINE VERIFIES 

C TOE SEQUENTIALINESS OF TOE DATA. 

C THE SUBROUTINE ALSO RETURNS HUMBER WHICH TELLS THE HIGHEST PLOT 

C HUMBER REQUESTED. THIS VALUE IS NEEDED FOR OTHER SUBROUTINES. 

C VARIABLE TOTAL CONTAINS TOE TOTAL NUMBER OF INDIVIDUAL PLOIS 

C BEING REQUESTED. 

C 

C THE VARIABLE K CONTAINS TOE NLSiBKR OF SEGMENTS THE DATA HILL 

C BE DIVIDED INTO. IT IS THE MAXIMUM NUMBER OF PLOTS PERMITTED. 

C 

C 

C PLOTS (100) CONTAINS I IF INDIVIDUAL PLOTB ARE DESIRED, 2 

C IF ALL SEGMENTS ARE TO BE PLOTTED. IN TOE LATTER CASE, 

C PLOTS HILL CONTAIN 1 TO K. 

C 

SUBROUTINE SELECT(PLOIS .NUMBER, K .TOTAL) 

INTEGER PLOTS (100), NUMBER, K.TOTAL.VER 


INITIALIZE 
50 NUHBER-0 

PRINT*, "COT TO SELECT ROUTINE" 

CHECK IF THE USER WISHES ALL POSSIBLE SEGMENTS OR JUST A FEW. 
PLOTS (100) IS 1 IF SELECTED PLOTTING IS DESIRED, 2 IF ALL ARE 
WANTED. 

IF (PLOTS (100) .HE. 1) COIO 400 

PRINT*, "ENTER TOE TOTAL NUMBER OF DESIRED PLOTS" 

READ*, TOTAL 

IF (TOTAL .EQ. K) GOTO 400 

PRINT*, "ENTER THE DESIRED PLOT NUMBERS, ONE AT A TIME, AND" 
PRINT*, "IN ORDER AS YOU ARE QUERIED." 

PUT THE PLOT NUMBERS INTO PLOTS . FIND THE HIGHEST PLOT 
NUMBER AND PUT IT IN NUMBER. 

DO 100 1-1 .TOTAL 

PRINT*, "PLOT NUMBER?" 

READ*,PLOTB(I) 

NUMBER-PLOTS ( I ) 

100 CONTINUE 

VERIFY THAT THE PLOT NUMBERS ARE lit ORDER. 

VER-0 

DO 200 1-1 .TOTAL 

IF (PLOTS (I) .LE. VER) GOTO 300 
VER- PLOTS ( I ) 

200 CONTINUE 

IF THE PROGRAM GETS HERE, VERIFICATION IS COHPLETE AND THE 
SUBROUTINE CAN RETURN TO WHERE IT WAS GALLED FROM. 

RETURN 

VERIFICATION ERROR COMES HERE. 

C 


300 PRINT*, «***PLOT NUMBERS ARE HOT SEQUENTIAL. TRY ACAD) « 
PRINT* , ‘|***RfHEHBER 1 IF YOU WANT ALL TIIR PLOTS ,REBP0ND w 
PRINT*, "TO T1IE NEXT QUESTION Win! ",K 
GOTO 50 


C 

• »{**»«« DE6IRED, PUT IN ALL POSSIBLE PLOT NUMBERS, 

C 

400 NUHBER-K 
TOTAL-K 
DO 500 I-I,K 
PLOTS (l)-I 
500 CONTINUE 
RETURN 
END 


c 

c 

c 


THE FOLLOWING SUBROUTINES WERE ADDED MARCH, 1984. 


SUBROUTINE MENU(QUES) 

INTEGER QUES 

THIS SUBROUTINE EXPLAINS THE VARIOUS PitOCE8BE8 WHICH CAN BE 
APPLIED TO STORED PLOTS. THE USER CHOOSES, SETTING QUES 
TO 1 “8 « 

PRINT*, "**********************iNnfrft#**ihMr##*#**'#***« i» 

PRINT*," " 

PRINT*," " 

SUBTR ACT PL0T #1 FROM ft. THE PLOTS MUST BE OF THE" 

prStvvork plot’" ™ E REEULTU!G PL0T 16 PLACED ® the “ 

PRINT*,’" " 

MAKE m INCREMENTED WINDOW AVERACE OP THE WORK " 
PRINT*’ ,,P "° T ' 71,8 aESU1,T 18 PLACED IN WORKPLOT." 

PRINT*, "3) TRADE TOE CONTENTS OF TWO PLOT BAY8. THE PLOTS" 
PRINT*, "NEED NOT BE OF THE SAME SIZE." 

PRINT*,"" 

AVEHACE PLOTS #1 AND ft TOGETHER. THE RESULT " 
PRINT*, "IS PLACED IN THE WORK PLOT." 

PRINT*," " 

PRINT*, "5) SHOW THE TIME, SIZE, AND PFT SIZE OF PLOT BAYS" 
PRINT*," " 

PRINT*, "6) PRINT OUT THE WORK PLOT ONTO THE SCREQI." 

PRINT*," " 

ra™«'i7ti G ° BACK P0!l MDHE ?L0TS °R QUIT THE PROGRAM.” 

rKIfiT" i 

PRINT*’" 8 " C0PY A ST0BAGE BAY 10 TOE WORK BAY." 

PRINT*,’" " 

PRINT* , "**»»*■»»»»»*»**»»* » * * ** ** * ********* * * » ♦**★«» 
PRINT*," " 

PRINT*, "TYPE IN YOUR CHOICE" 

READ*, QUES 
PRINT*," " 

PRINT*," " 

RETURN 

END 


BUBROUTINE DIPFSET<SAVMAOl , BAVFRQl , P ARAMS 1 .8AVMAG2 >8AVFRQ2 . 

* PARAN52 .WORKHAC, WORKFRQ, WKPARMS ) 

THIS SUBROUTINE SUBTRACTS PI,OT#l FROM PLOTft. THE TWO PLOTS MUST 
BE OF THE SAKE 6IZE, THE RESULTING SPECTRUM IS PUT IN WORKKAG. 

INTEGER qtfES 

REAL JAVMaGI ( 1025) , SAVHAC2 (1 025 ) , SAVFRql (1025 ) , 

* BAVFRQ2 (1025 ) , WQRKMAC( 1025) ,WORKPRq( 1025), 

* PARAK31 (10) ,PARAMS2(10) ,WKPARKS(10) 

PRINT*,"THIB ROUTINE SUBTRACTS PLOT01 FROM PLOT02, T1IE PLOTS" 
PRINT*, "MU8T BE OF EqUAL SIZE. THE RECULTDiC PLOT IS" 

PRINT*, "PLACED IN THE WORK PLOT" 

CALL PRTSTAT(PARAKS1 , PARAMS 2 , WKPARMS) 

PRINT*, "DO YOU NEED TO SWITCH PLOTS 1 AND 2? 1-YES" 

PRINT*," 2 -NO" 

PRINT*," 3-ABORT ROUTINE" 

READ* ,QUES 

IP (qUES ,Eq. 3) RETURN 

IF (QUES .EQ. 1) CALL TRNSPLT(SAVMACl ,SAVFRQl ,PA£lAH81 .8AVMAC2 , 

* 8AVFRQ2 ,PARAMS2 ) 

CALL DIFFRNC(SAVMAC1 .BAVFRQl , PAR AMS 1 .8AVMAC2 .8AVFRQ2 .PARAKS2 , 

* WORKKAG, WORKFRQ ,WKPARM8 ) 

RETURN 

END 


SUBROUTINE WND08E T ( SAVMAG1 , BAVFRQl .PARAM31 .8AVMAC2 .8AVPRQ2 , 

* PARAMS2 , WORKMAC , 

* WORKFRQ , WKPAR MB ) 


REAL BAVHAC1 (1025) .SAVMAC2 (1025), SAVFRQ1 (1025) ,SAVFRQ2( 1025 ) , 

* PARAMSI (10},PARAHS2(10) ,WORKHAC(1025) ,WQRKFRq(l025) , 

* WKPARHS(IO) 

INTEGER QUES 


THE FOLLOWING DOES AN INCREMENTED WINDOW AVERAGE OF THE WORK PLOT. 
THE RESULT IS PUT IN WORK PLOT ALSO, 


PRI>tf'*,"THIS ROUTINE COMPUTES AN INCRQ1ENTED WINDOW AVERAGE OP THE ” 
PRINT*, "WORK PLOT. THE RESULT IS ALOS PUT INTO THE WORK PLOT." 
PRINT*, "DO YOU NEED TO PUT A SAVED PLOT INTO THE WORK PLOT," 

PRINT*," 1-YEG,DAY01" 

PRINT*, " 2-YES , BAY 02" 

PRINT*," 4-NO" 

PRINT*," 5-ABORT" 

READ*, QUES 

IF(QUES .EQ. 1) CALL TRNSPLT(SAVKAG1 , BAVFRQl , PARAMSI , WORKKAG, 

* UORKFKQ , WKPARMS ) 

IF(QUES . EQ. 2) CALL TBNSPLT(SAVMAG2 .SAVFRQ2 , PAR AMS 2 .WORKKAG, 

* WORKFRQ, WKPARMS) 

IF(QUES .EQ. 5) RETURN 

CALL WNDQAVE(WORKMAG, WORKFRQ, WKPARMS) 

RETURN 

END 


BUBROUTINE T8NSBET(8AVMAC1 , BAVFRQl .PARAMSI .8AVMAG2.SAVFRQ2, 

* PARAM82 , WORKKAG, 

* WORKFRQ, WKPARMS ) 

REAL SAVMAG1 (1025 ) , 8AVMAG2 ( 1025 ) , BAVFRQl ( 1025 ) , SAVFRQ2 ( 1025 ) , 

* WORKMACd 025 ),W0RKFRQ( 1025) .PARAMSI ( 10) ,PARAMS2( 10) , 

* WKPARMS (10) 

INTEGER QUES 
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c switch plots around 

c 

PRINT*, "THIS ROUTINE TRADES PLOTS BETWEQJ A SPECIFIED PLOT" 
PRINT*, "BAY AND THE WORK PLOT." 

PRINT*, "W1IICU PLOT DO YOU WISH TO TRADE WITH THE WORK PLOT" 
PRINT*," < 4-NO SWITCHING )" 

READ*, QUES 

IF (qUES .EQ. 1) CALL TRN6PLT(8AVMACl , SAVFRQl ,PARAMSl ,VORKHAG, 

* WORKPRQ , WKP ARMS ) 

IF (OUE9 .EQ. 2) CALL TRNSPLT(SAVMAC2 .SAVFRQ2 .PARA1IB2 .WORKHAG, 

• WORKFRQ , WKPARHS ) 

CALL PHTSTAT(rAPAHSI , PAR AMS 2 .WKPARHS ) 

RETURN 

END 


SUBROUTINE C0PY8ET(SAVHACl .SAVFRQl .PARAHSl ,SAVHAC2 ,SAVFRQ2 , 

* PARAHS2, WORKHAG, WORKFRQ, WKPARHS) 

C THIS ROUTINE COPIES A STORAGE BAY TO THE WORK BAY. NOTE THAT 

C THIS WILL DESTROY WHATEVER IS IN THE WORK BAY PREVIOUS TO THE 

C CALL. TO COPY THE WORK BAY TO A STORACE BAY, TRADE THE WORK 
C BAY WITH THE DESIRED STORACE BAY, AND THEN DO THE COPYING. 

C 

INTEGER QUES 

REAL SAVMACl ( 1025 ) , BAVHAC2 Cl 025 > , WORKMAGC1025) , 

* SAVFRQ2( 1025) ,WoRKFRQ(I 025) , PARAHSICIO) , PARAHB2C 10) , 

* WKPARHS (10) 

PRINT*, "COPIES A STORAGE BAY TO THE WORK BAY" 

PRINT*," " 

PRINT*, "WHICH 
PRINT*," 

PRINT*," 

PRINT*/* 

READ*, QUES 
IF (QUE3 .EQ. 

* 

IF (QUES .EQ. 

* 

IF (QUES .EQ. 

RETURN 
END 


SUBROUTINE AVRGBETCSAVHACI .SAVFRQl , PAR AMS 1 ,SAVHAG2,SAVFRQ2, 

* PARAMS2 .WORKHAG, WORKFRQ, WKPARHS) 

REAL W0RKMAG( 1 025 ) ,WDRKFRQ( 10J5) ,SAVHAG1( 1025) ,SAVMAG2(10'5) , 

* SAVFRQl ( 102 5) ,SAVFRQ2( 1025) 

IN TEGER WKP ARMS ( 1 0 > , PAR AMS 1 ( 1 0 ) , PAR AMS 2 ( 1 0 ) , QUES , ERR 
PRINT*, "AVERAGES PLOT BAYS I AND 2. 1-CONTINUK" 

PRINT*," 2-ABORT ROUTINE" 

HEAD*, QUES 

IF (QUES .EQ. 2) RETURN 
C 

C MAKE SURE BAYS 1 AND 2 ARE OF THE SAME SIZE 
C 

ERR-0 

CALL SIZEERR(PARAHSi,PAHAKS2,2,ERR) 

IF(ERR .NE. 1) RETURN 
C 

C PERFORM AVERAGING 
C 

CALL AVERACE(BAVHAGl .SAVFRQl ,PARAMSl .SAVHAG2 .SAVFRQ2, 

* PARAMS2.WORKMAG, WORKFRQ, WKP ARMS) 

RETURN 

END 


STORAGE BAY DO YOU WISH TO COPY?" 

1- bay tfi" 

2- BAY 02*' 

3- ABORT ROUTINE" 

1) CALL UOPY(SAVMAGl, SAVFRQl, PARAMSl .WORKMAC, 

WORKFRQ, WKP ARMS) 

2) CALL COPY(SAVMAG2,SAVFRQ2,PARAHS2 J WOR»UG, 

WORKFRQ, WKP ARMS) 

3) RETURN 
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SUBROUTINE DIFPRNCCLOCHAGl ,L0CFRQl , LOCP/Jtl , L0CHAC2 , L0CPRQ2 , 
* L0CPAR2 .WORKMAC , WORKFRQ , WKPARMS ) 


INTEGER I,ERR,K 

REAL LOCMAG1 (1025), IOCFRQ1 (1025), L0CHAG2<1025) ,LOCPRQ2(1025), 

* WORKMAG(1025) ,W0RKFRQ(1025),L0CPAR1(10) ,LOCFAR2(10) , 

* WKPARMS (10) 

CHECK TO SEE THAT THE TWO ARRAYS HAVE SAKE SIZE. 

ERR-0 

CALL 8IZEERR(LOCPARl,LOCPAR2,l,ERR> 

IP (ERR .RE. I) RETURN 

TAKE THE DIFPERH1CE 

K“IFIX(L0CPAR1 (2) + .1) 

DO 100 Z-1,K 

WORKMAG( I ) -LOCHAC2 ( I )-LOCHACl ( I ) 

W0RKFRQ(I)-L0CFRQ1(I) 

100 CONTINUE 

PUT LOCPAR1 INTO WKPARMS, I.E. EXCHANCE PARAMETERS 

WKPA3HB(l)-0.0 
WKPARKS ( 2 )=L0CPAR1 (2 ) 

WKPARMS (3 ) -LOCPAR1 ( 3 ) 

RETURN 

EDO 


SUBROUTINE WNIIOAVE ( WORKMAG , WOfJKFRQ , WKPARHS ) 

INTEGER TOTLFRQ ,K1 ,N , BINS , EN DCNT 

REAL WNDOS IZ , POWRLAW , W DOFRQ( 1 02 5 ) , Fl , HLDAB BY ( 1 025 ) , 

* WORKKAG ( 1 02 5 ) , W0RKFRQ( 1 025 > , MULTS IZ , WKPARKS (10) 

PRINT*, "THIS ROUTINE COMPUTES A WINDOWED AVERAGE OF A" 

PRINT*. "PLOT IN THE WORK PLOT." 

PRINT*, "THE RESULTING PLOT IS PUT IN THE WORK PLOT." 

PRINT*, "IT WILL HAVE A DIFFERENT SIZE THAN THE INPUT" 

PRINT*, "PLOT." 

PRINT*, "THE FIRST FREQUENCY OF THE INPUT PLOT IS ",WORKFRQ{2) 
PRINT*, "WHAT SHOULD THE FIRST WINDOW FREQUENCY BE?" 

WNDOFRQ(l) SHOULD BE 0.0 , CORRESPONDING TO WORKFRQ'S FIRST 
FREQUENCY BEING 0 .0 . 

READ*,WNDOFRQ(2) 

WNDOPRQ(I)-0.0 
PRINT*," " 

VNDOSIZ-WNDOFRQ(2) 

PRINT*," " 

PRINT*, "HOW MANY WINDOWS DO YOU WANT?" 

READ* , TOTLFRQ 
PRINT*," " 

PRINT*, "WHAT SHOULD THE MULTIPLICATIVE WINDOW INCREMENT EE?" 
PRINT*, "(SHOULD BE GREATER THAN 1)" 

READ*, MULTS IZ 
PRINT*," " 

PRINT*. "TO WHAT POWER LAW DO YOU WISH TO FIT THE SPECTRUM?" 
PRINT ‘."(SHOULD BE A NEGATIVE NUMBER)" 

READ*, POWRLAW 
PRINT*," " 

VALIDATE POWRLAW AND MULI5IZ 

IF (POWRLAW .GT. 0.0) CALL PRNTERR(I) 

IF (POWRLAW .GT. 0.0) RETURN 
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IF (HULISIZ .LT. 1,0) CALL PRHTERR(2) 

IF (HULISIZ ,LT. 1.6) RETURN 

SET UP AN ARRAY OF THE LOWER FREQUENCIES OF EACH WIDOW 
THE FIRST TWO FREQUENCIES WERE DETERMINED BY WHDOSIZ AND " 
THE USER SPECIFIED FIRST FREQUENCY, 


DO 50 I«3,?f)TLFRQ 

WHDOFRQ( I MWND0FRQ( I"1 )-KHD0FRQ( I-2))*HULT5 IZ»WNDOFRQ( 1-1 ) 

50 CONTINUE 

HOW SEARCH THROUGH THE WORKFRQ, FIND ALL POINTS THAT FALL IN 
A BIN, AVERAGE THEM, AND PLACE EACH AVERAGE SUCCESSIVELY IN 
THE WORKHAG. IF NO FREQUENCIES FALL IN A BIN, PUT ZERO IN THE 
PROPER WORKHAG SLOT. 

THIS HICHT PLACE ZERO IN THE FIRST BIN AND IN THE UPPER BINS 
IF TOTLFRQ TAKES THE SEARCH ABOVE THE HIGHEST FREQUENCY IN 
THE WORKFRQ ARRAY. 

THERE IS ONE LESS BIN THAN THERE ARE WINDOW FREQUENCIES. 

SO BINS“T0TLFRQ-1 . 

BIN8-T0TLFRQ-1 

LOOK THROUGH ALL BINS FOR APPROPRIATE POINTS 

ALSO WATCH FOR Q1PTV BINS. AFTER 2 Q1PIY BINS IT CAN BE ASSUMED 
THAT NO PREQUECIES EXIST IN THE HIGHER BINS. COUNT THESE WITH ENDCNT. 
EVERY TIME WE GET A BIN WITH AT LEAST ONE COUNT, RESET THE 
ENDCNT TO 0 

K1«IFIX(WKPARHS(2) .1) 

ENDCHT-0 
DO 300 I-l.BINS 
N-0 

ADDUP-0 .0 
DO 200 J«1,K1 

WNDOFRQ(I) IS THE LOWER FREQUENCY OF THE I'TH BIN, AND 
WNDOFRQ(Ul) IS THE UPPER FREQUENCY OF THE I'TH BIN 

IF((WORKFRQ(J) .CT. WNDOFRQ(I)) .AND. (WORKFRQ(J) .lt. 

* WNDOFRQ(m))) 110,200 

110 N-N+l 

ADDUP«ADDUP*WORK"!,;C(J) 

200 CONTINUE 

AT MlIS POINT, N»0 MEANS NO POINTS WERE FOUND. 

ADJUST ENDCNT ACCORDINGLY. 

IF(N .EQ. 0) HLDARRY(I)-0.0 

IF(N .HE. 0) ENDCNT*>0 

IF(N .EQ. 0) ENDCNT-ENDCNT+I 

IF (ENDCNT .EQ. 2) WKPARHS(2)»FL0AT(I) 

IF<N .NE. 0) HLOARRY(I)«ADDUP/ (FLOAT(N) ) 

300 CONTINUE 

CHECK IF THE DATA BINS WERE NOT ENOUGH TO COVER THE INPUT 
DATA. DO THIS BY CHECKING FOR ENDCNT LESS THAN 2. IF THIS 
IS THE CASE THEN WKPARHS(2) SHOULD BE SET TO BINS. 

IN ANY CASE PRINT OUT HOW MANY BINS WERE ACTUALLY USED. 

IF (ENDCNT .LT. 2) WKPARMS(2)»FLOAT(BINS) 

PRINT*, "TH ERL ARE ",WKPARMS(2) FREQUENCY BINS IN THE" 

PRINT*, "INCRaiENTED WINDOW AVERAGE SPECTRUM." 

NOW WE HAVE A COLLECTION OF AVERAGE MAGNITUDES FOR A SET OF BINS. 

TO PROPERLY ASSIGN EACH VALUE TO A PARTICULAR FREQUENCY IN 
A BIN WE MUST ACCOUNT FOR THE POWER LAW WEIGHTING. 

DO 400 I“2 ,BINS 

FI “WNDOFRQ(I) -WND0FRQ( 1-1 ) 
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WORKFRQ<I)-(UWHDOFRQ(l)*<MULTSlZ+l,0)-MULTSIZ*VNDOFRQ(I-l))** 

* (POWRLAWfl .0)-WND0FRQ(l)**(P0WRtAW+l .0) )/ 

* ( (POWRLAWtl ,0)*MULTSIZ*(Fl) ) )**(1 .0/POWRLAW) 

400 CONTINUE 


NOW PUT THE TEMPORARY HOLDING ARRAY (HLDARRY) INTO THE 
WORKMAG ARRAY. 

El“-mX(WKPARHS(2) 4- .1) 

DO 350 1-1 ,KI 
WORKMAG (I } »HLDARRY< I ) 

350 CONTINUE 

CHANGE THE PARAMETERS IN WKPARMS 
WKPARMS(3) THE FFT SIZE IS THE SAME. 

WKPARMS (1) -0.0 

RETURN 

END 


T11E FOLLOWING SUBROUTINE MOVES ONE WORK BAY TC ANOTHER 


SUBROUTINE TRNSPLT(L0CMAG1 .LOCFRQl .L0CPAR1 .LOCHAG2 ,L0CFRQ2 , 

* L0CPAR2) 

REAL LOCMAG1 ( 1 025 ), LOCMAG2( 1025 ),U>CFRQl( 1025), 

A L0CFRQ2( 1025 ) , L0CPAR1 ( 1025 ) , LOCPAR2 ( 1025 ) , 

* SAVEMAG,SAVEFSq,SAVEPAR 

INTEGER I,J,K,HAXSIZE 

FIND WHICH OF THE TWO ARRAYS IS LARGER AND USE THE INTEGER 
VALUE IN THE LOOP. 

HAXSIZE«IFIX(LOCPARl(2) + .1) 

IF (LOCPAR2(2) ,GT. LOCPARl(2)> KAXSIZE“IFIX(L0CPAR2(2)+.l) 
SWITCH ELEMENTS OF THE TWO ARRAY SETS 

DO 200 I-i.MAXSIZE 
SAVEHAC-LOCMAGI(I) 

LOCMAGl ( I ) -LOCMAC2 ( I ) 

L0CHAC2(I)-SAVEMAG 
SAVEFRQ-L0CFRQ1 ( I ) 

LOOFRQI ( I ) !B L0CFRQ2(I ) 

L0CFRQ2 ( I ) -SAVE FRQ 
200 CONTINUE 

DO 300 1-1,10 
SAVEPAR-LOCPARl(l) 

LOCPAR 1 ( I ) -L0CPAR2 U ) 

LOC PAR 2(1) -8 AVEPAR 
300 CONTINUE 

RETURN 

END 

SUBROUTINE COPYCLOCKACl ,L0CFRQ1 .IOCPARl .L0CMAC2 ,L0CFRQ2, 

* L0CPAR2) 


THIS ROUTINE COPIES LOCAL BAY#1 TO LOCAL BAY #2 . THE CONTENTS 
OF BAY #2 ARE DESTROYED. 

INTEGER I, SIZE 

REAL L0CHAG1( 1025 ) , LOCHAC2 0025) j LOCFRQl (1025 ) , 

* LOCFRQ2(1025),LOCPAR1(10),LOCPAR2(IO) 


S1ZE“IFIX(L0CPAR1(2) + .1) 
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DO 200 1-1, SIZE 

LOOHAC 2 ( T. ) -LOCHAC 1 ( I ) 

LOCFRQ2(l) -LOCFRQl ( I) 

200 CONTINUE 

DO 300 1-1,10 

LOCPAR 2 ( I ) -LO0PAR 1 U ) 

300 CONTINUE 

RETURN 

END 

C 

C SUBROUTINE FRTSTAT GIVES TOE PLOT BATS AND TltEIR SIZES 
C 

SUBROUTINE PRTS TAT (PARAHS 1 , PAR AMS 2 .WKPARMS) 

REAL PAR AKB 1 ( 10 ) , PARAMS 2(10), WKPARMS (10) 

INTEGER I , J , K 

PRINT*, "PLOT BAT THE SIZE (POINTS) FFT SIZE" 
PRINT* , "************************«********************" 
PRINT*," " 

WRITE(2,900) l,PARAMBl(i) ,PARAMS1(2),PARAMS1(3) 

WRITE(2 ,900) 2 ,PARAMS2(1) ,PARAMS2(2) ,PARAMS2(3) 

WRITE(2 ,1000) WKPARMS Cl ) .WKPARMS (2) .WKPARMS (3) 

900 FORMAT(I1,11X,F6.3,2X,F5.0,13X,F5.0) 

1000 FORMAT(9HWORK PLOT,3X,F6.3,2X,F5.0,13X,F5.0) 

RETURN 

END 


C 

C SUBROUTINE, AVERAGE TAKES THE AVERAGE OF TWO PLOTS . 

C 

SUBROUTINE AVERACE(L0CMAG1 .LOCFRQl .LOCPARl .L0CMAG2 , 

* LOCFRQ2 , L0CPAR2 .WORKMAG, WORKFRQ .WKPARMS ) 

REAL L0CMAG1 ( 1025 ) , LOCMAG2 ( 1025 ) ,W0RKMAG( 1025 ) , LOCFRQl ( 1025 ) , 

* LOCFRQ2 ( 1025 ) , WORKFRQ( 1025 ) , LOCPARl (10) .LOCPAR2 ( 10) , 

* WKPARMS(IO) ,ADDUP 
INTEGER I, J, POINTS 

POINTS“IFIX( LOCPAR 1(2) + .1) 

DO 100 1-1 .POINTS 

W0RKHAG( I ) - (LOCHAG1 ( X)+LOCHAG2 (I) ) /2 .0 
W0RKFRQ( I > “LOCFRQU I ) 

200 CONTINUE 
C 

C PUT PARAMETERS INTO TOE WORK PARAMETERS 
C 

WKPARMS CD— 0 ,0 
WKPARMS ( 2 ) -LOCPAR 1 ( 2) 

WKPARMS (3 )-LOCPARl (3) 

RETURN 

END 

C 

C THE FOLLOWING TWO ROUTINES DEAL WITH ERRORS 
C 

SUBROUTINE PRNTERR(l) 

INTEGER I 
<7\T0(10,20,30)I 
10 PRINT*," " 

PRINT*,"*** POWER LAW MUST BE A NEGATIVE NUMBER***" 

PRINT*," " 

RETURN 

20 PRINT*," " 

PRINT*, "***THE MULTIPLICATIVE INCREMENT OF THE WINDOW MUST" 
PRINT*," BE GREATER THAN 1.00***" 

PRINT*," " 

RETURN 
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30 PRINT*, " " 

PRINT*, "***THERE ARE ONLY EIGHT ALLOWED COMMANDS***" 
PRINT*," " 

RETURN 

END 


SUBROUTINE S1ZEERR{L0CPARI .LOCPAR2 ,I,ERR) 


REAL LOCPARl(IO) ,LOCPAR2{IO) 
INTEGER I , ERR , SIZE1 ,SIZE2 


THIS ROUTINE CHECKS IF TWO ARRAYS HAVE THE SAME SIZE. ERR-1 
IF THE ARRAYS HAVE THE SAME SIZE. OTHERWISE ERR-0 . 


SIZEl=IKIX(LOCFARl(2) + .1) 

SIZE2-IFIX(L0CPAR2(2) + .1) 

IF (SIZE1 ,EQ. SIZE2) ERR-1 
IF (ERR .EQ. 1) RETURN 
GOTOClO ,20)1 
10 PRINT*," « 

PRINT*, "***ARRAYS ARE OF DIFFERENT SIZES IN DIFFRNC***" 
PRINT*," " 

RETURN 

20 PRINT*," " 

PP.INT*, "*** ARRAYS ARE OF DIFFERENT SIZES IN AVRGSET***" 
PRINT*," " 

RETURN 

END 
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11.5 DATPLOT 


PROGRAM DATPLOT (INFUT, OUTPUT, DATA29 ,TAPE1 •INPUT, 

♦TAPE2-OUTPUT ,TAPE3 “MTA29 > 

C 

c THIS PROGRAM READS 2000 POINTS FROM A ROCKET FILE AND PLOTS 
C THE POINTS AS A FUNCTION OF TIKE. THIS PARTICULAR VERSION ASSUMES 
C A RECORD LEKCTH OF 2 SO (SEE OTHER ROCKET PROGRAMS FOR ROCKET FORMAT) 
C AND A DIGITIZING RATE OF 5000. THE PLOTTING STARTS AT STRT1M 
C SECONDS INTO THE FILE, RELATIVE TO THE FILE'S BECINtUNC. 

C THIS ALLOWS USE WITHOUT KNOWLEDGE OF THE ABSOLUTE TIME OF EACH 
C RECORD. THE USER MUST ONLY KNOW THE DESIRED POSITION RELATIVE 
C TO THE PARTICULAR FILE'S STARTING TIME. 

C TO RUN THIS PROGRAM THE CCS PLOTTING LIBRARY MUST BE ACCESSED 
C WITH! CRAB,GC8(TEKT,3D). THIS IS ENTERED PRIOR TO RUN. THE 
C LIBRARY IS NOT NEEDED TO COMPILE DATPLOT. 

C 

REAL AR RAY (1251), TIME ( 2000 ) , CHAN ( 2000 ), TSTPRNT , STRTIH 
INTECER IARRAY(1250) ,1 , J,K, L, BTRTREC , DICRATE , SKIP , MAXMAG , 

* QUES 

EQUIVALENCE (ARRAY(2) ,XARRAY(1) ) 

C DIGRATE IS THE DIGITIZING RATE OF THE DATA. 

DIGRATE** 5000 

C MAXMAG IS THE LARGEST EXPECTED VALUE OF DATA. 

MAXMAG** 256 

FRINT*, "THIS PROGRAM PLOTS 2000 POINT8. IT ASSUMES" 

PRINT*, "A ROCKET FILE AS THE INPUT." 

PRINT*," " 

PRINT*, "AT HOW MANY SECONDS INTO THE PILE SHOULD THE PLOT BEGIN?" 
READ*, STRTIH 
CALL USTART 
30 REWIND 3 
C 

C STRTREC IS THE FIRST RECORD WHICH WILL BE PLOTTED. 

C THE FOLLOWING LOOP PUTS THE FILE READ POINTER AT THE PROPER PLACE. 

C THE RECORD IS NOT READ. 

C THE "+ .1" HELPS PREVQIT ROUNDOFF ERRORS FROM THE IFIX 
C 

STRTE EC"IFIX( (FLOAT (DIGRATE/250 ) )*STRT1M t .l)+i 
SKIP-STRTREC-1 
IF (SKIP .EQ. 0 ) GOTO 50 
DO 40 L-l.SKIP 
READ (3) ARRAY 
40 CONTINUE 
50 DO 300 K“l,8 
READ (3) ARRAY 
TIME ((K-1)*250+1)“ARRAY(1) 

DO 100 1-2,250 

TIME ( I* (K-l)*250) “TIME (1+ (K-l ) *250-1 )+(l ,0/ (FLOAT(DIGR.IE) ) ) 

100 CONTINUE 
C 

C NOW LOAD UP THE ARRAY CALLED CHAN. NOTE THAT CHAN IS REAL AS 
C ALL GCS GRAPHICS ROUTINES REQUIRE REAL NUMBERS AS INPUTS. 

DO 200 J-1,250 

CHAN ( J+ (K-l)*250 )»FLOAT(lARRAY( J*5-3 ) ) 

200 CONTINUE 
300 CONTINUE 
C 

C SET UP A USER DEFINED AREA THE SIZE OF A TEKTRONIX SCREBJ. 

C ALSO SET THE VIRTUAL UNITS IN THE RANGE OF THE DATA BY MEANS OF A 
C UWXNDO COMMAND. 

C 

CALL UPSET ( "PRECISION" ,-5 . ) 

CALL USETCOWNSCALE") 

CALL UDAREA(0 .0,7. 0,0. 0,5. 7) 

CALL UWIND0(TIME(1 ) ,TIME(2000) ,0 .00, FLOAT (MAXMAG) ) 

CALL UELOTK TIME, CHAN, 2000.0) 

CALL UPAUSE 

PRINT*, "DO THE NEXT ONE? 1-YES" 

PRINT*,"' 2-NO" 
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READ*,QUE8 

CALL UERASE 

STRTIH-BTRTIM+.4 

IP (qUES .EQ. 1 )C0T0 30 

CALL UEND 

STOP 

END 
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II. 6 SPIN FIX 


PROGRAM SPINFIX(INPUT,0UTPUT,DATAIN,DAT0UT,TAPE1"IKPUT, 
♦ TAPE 2 ■OUTPUT , TAPE 3 “DATA IN .TAPE4-DATOUT) 


PROGRAM BPINPIX BY RODNEY STOLTZPUS 11/1/83 


THIS PROCHAH PILTERS OUT THE SPIN FREQUENCY 
FROM ROCKET DATA. THIS FREQUENCY IS OF THE 
ORDER OF 4-8 HZ. THIS FILTERING INVOLVES CONVOLVING THE 
THE DATA WITH A RECTANGULAR WINDOW 1 HIGH AND OF LENGTH 
C SPIN PRD. THIS CONVOLUTION TRANSFORMS TO MULTIPLYING THE 

C SPECTRUM DY A SIN SQUARED FUNCTION, PUTTING HOLE IN THE 

C SPECTRUM EVERY 1/BPINPRD HZ. THE CONVOLUTION PERFORMED IS 
C NOT EXACTLY CONVBmONAL. INSTEAD OF ADDING UP THE WHOLE 

C INTERSECTING GROUP OF POINTS , THE LEAD POINT IS ADDED AND THE 

C TRAILING POINT IS SUBTRACTED FROM THE PREVIOUS SUM. THE 
C INITIAL VALUE IS SET AT ZERO. THE ONLY EFFECT THIS HAS IS TO 

C DISTORT THE OFFSET. ALSO, EACH POINT IS NOT NORMALIZED BY 

C THE TIME INCREMENT BETWEEN POINTS. 

C 

C ARRCNT SUBBLOCKS OF DATA ARE READ INTO A LARCE ARRAY. 

C THIS DATA IS ANALYZED UNTIL ANOTHER SUBBLOCK IS NEEDED, 

C WHEREUPON DATA IS READ IN AND PLACED IN THE APPROPRIATE 
C POSITION WITH RESPECT TO THE OTHERS. THUS THE 
C DATA CYCLES THROUGH THE MAIN BLOCK. 

C THE MAIN BLOCK MAY BE VIEWED AS CIRCULAR. TWO WIPERS 

C ROTATE, SAMLING THE DATA- WIEN THE TRAILING WIPER 

C FINISHES A SIIBBLOCK A NEW ONE IS INSERTED. 

C THE COUNTER REDONT RUNS FROM 1 TO HOWEVER MANY DATA 

C POINTS ARE DESIRED. THIS WILL RUN INTO THE THOUSANDS. 

C EVERYTIHE IT PASSES AN INCREMENT OF ARRCNT*250 THE FUG OVURFLW 
C IS INCREMENTED. THE COUNTER PNTCNT INDEXES THE SUBBLOCK 
C IDATA. WHEN IDATA IS FULL IT IS WRITTEN INTO TAPES. 

C PNTCNT IS RESET AFTER A DATA RECORD IS WRITTEN OUT. . 

C 

C WHEN RUNNING THE OUTPUT DATA THROUGH THE "SPECTRA "PROGRAM 
C YOU MUST MULTIPLY THE SPECTRUM BY (CONSTANT )*{'FRE- 

C QUENCY)**2 TO ADJUST FOR THE RECTANGULAR WINDOW. 

C 

C 

C IMPORTANT NOTE 

C THIS PROGRAM ASSUMES DATA RECORDS IDENTICAL TO THOSE OF 

C ROCKET SHOTS 31.029 AND 31.028. THAT MEANS A DIGITIZING 

C RATE OF 5000, RECORD SIZE OF 1250, AND FIVE DATA CHANNELS 

C PER FRAME , THE SECOND CHANNEL BEING THE ONE OF INTEREST 

C 

C VARIABLE DICTIONARY 

C 

C TIMEL1 THE TIHE MARKER OF THE RECORD PRIOR TO THE 

C STJRTING RECORD. 

C 

C BLKCNT THE NEXT SUBBLOCK TO BE FILLED. 

C 

C ARRAY CONTAINS 1251 DATA POINTS. AS USUAL THE 

c FIRST CONTAINS THE CORRESPONDING TIME. 

C 

C IARRAY IS THE EQUIVALENT OF ARRAY BUT IN INTEGER 

C FORK . 

C 

C IRAWDAT CONTAINS FOUR BUBBLOCKS WITH THE TIMF. AND THE 

C FOUR OTHER DATA CHANNELS STRIPPED OFF. 

C 

C SPINFRQ IS THE SPIN FREQUENCY OF THE ROCKET. 

C 

C 8PNPRD IS THE PERIOD OF THE SPIN IN TERMS OF 

C DATA POINTS. 

C 

c ARRCNT is THE NUMBER OF SUBBLOCKS IN THE MAIN BLOCK. 
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C THEBE SHOULD BE EMOUGH SUBBLOCKS TO COVEB THE SPIN PERIOD 

C PLUS AN EXTRA SUB BLOCK. 

C 

C LEFTOVR IS THE DIFFERS) CE BETWEEN THE BLOCK LQiGTH AND 

C THE SPIN PERIOD. 

C 

C IDATA CONTAINS THE INTEGER FORM OF THE PROCESSED DATA. 

C 

C I, I ARE ASSORTED COUNTERS. 

C 

C NORM IS THE NORMALIZED INDEX (POSITION IN IRAWBAT) 

C OP THE DATA POINT SPNPRD AHEAD OF THE CURRENT POINT 

C OVERFLW TELLS THE INCREMENT OF ARRCNT*250 THAT REDCNT 

C IS OK. SEE PROGRAM NOTES ABOVE THIS DICTIONARY. 

C 

C REDCNT INDEXES THE ARRAY IRAWDAT. REDCNT IS UNBOUNDED. 

C SEE OVERFLH. REDCNT COUNTS THE POINTS THAT HAVE 

C BEEN READ ALREADY. 

C 

C PNTCNT INDEXES THE ARRAY IDATA. PNTCNT IS CYCLED FROM 

C 2 TO 251 . 

C 

c 

REAL ARRAY (I 251 ) ,SPINFRQ,STRTIHE , TIME, TIMEL1 .SPAN 
INTEGER SPNPRD , ARRCNT , LEFTOVR , IDATA(251 ) , IRAWDAT( £000 ) , 

+ I , J, BLKCNT ,REDPNT , PNTCNT .OVERFLH , IARRAY (1250) , SPANPTS 

EQUIVALENCE (ARRAY(2 ), IARRAY( L) ) 

INPUT THE SPIN FREQUENCY. 

PRINT*, "WHAT IS THE SPIN FREQUENCY OF THE ROCKET" 

READ*, S PIN FRQ 

PRINT*, "WHAT IS THE STARTING TIME IN U.T. SECONDS" 

READ*, STRTIHE 

PUT THE DATA CORRESPONDING TO STRTIHE INTO ARRAY. 

10 READ(3> ARRAY 

TIMEL1-STRTIME-0 .05100. 

IF (ARRAY(l) ,GT. TIHEL1) GOTO 20 
GOTO 10 

20 CONTINUE 

TIME-ARRAY(l) 

FIND BPIN PERIOD, BLOCK SIZE, AND LEFTOVER. 

8PNPRD-IFIX(5000 .0/SPIN FRQ) 

ARRCNT"2+((IPIX( (5000 .0/SPINFKQ)*. 5)1/250) 

LEFTOVR- (ARRCNT*250)-SPNPRD 

LOAD UP IRAWDAT WITH ROCKET DATA 

DO 100 1*1,250 
IRAWDAT ( I ) -IARRAY ( ( 1-1 ) *5+2 ) 

100 CONTINUE 

DO 120 J«2,ARRCirr 
RRAD(3) ARRAY 
DO 110 1-1,250 

IRAWDAT( 1+ ( J- 1 ) *2 50 ) -IARRAY ( ( 1-1 ) *5+2 ) 

110 CONTINUE 
120 CONTINUE 

BEGIN FILTERING 

BLKCNT-1 

IDATA(1)»0 

REDPNT-1 

PNTCNT-2 

OVERFLW-O 

ASK FOR THE DESIRED NUMBER OF SECONDS. 

MUST BE A MULTIPLE OF 0.05000 . 
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PRINT*, "HOW MANY SECONDS OK DATA DO YOU WANT?" 

READ*, SPAN 

SPANPTB-IFIX(8PAN*5000 .00*. 1 ) 

DO 300 I-1,8PANPTB 

CALL PROCESS (PNTCNT .REDPHT , IRAWDAT, IDATA .SPKPRD.TIKE , 
♦ OVERFLV,BLKCHT,ARRCHT) 

300 CONTINUE 
STOP 

mo 


SUBROUTINES CETNEXT AND PROCESS 


SUBROUTINE CElNEXT(BLKCNT,LOCIRAW) 

READ LOCARRY(125l) 

INTEGER LOCIARR(1250) ,1 ,J,LOCIRAW(2000) ,BLKCNT 
EQUIVAI.aiCE(L0CARRY(2} , LOCIARR(l) ) 

READ(3) LOCARRY 

PUT DATA INTO L0CIRAW AT THE POSITION BLKCNT. 
DO 400 1-1,250 

LOCIRAW( ( BLKCNT-X ) *2 50 +1 ) -LOCI ARR( (I- 1 ) *5*2 ) 
400 CONTINUE 
RETURN 
END 


SU BROUT IKE PROOES S ( LPNTCNT .LREDPNT , IRAWDAT , IDATA , 

+ S PNPRD, TIME, OVER FLW, BLKCNT ,ARRCNT) 

REAL OUTPUT(1251),TIME 

INTEGER IOUTPUTd 250 ) , LPNTCNT .LREDPNT , IRAVDAT(2000) , IDATA(251 ) , 
♦ 8PNPRD, OVERFLW , BL KC NT , NORM , OFFSET , AR RCNT 

EQUIVALENCE (OUTPUT< 2) , IOUTPUT(l ) ) 

IF ((LREDPNT) (ARRCNT*250*(0VERFLW+1))) OVERFLW- 
+ OVERFLH* 1 


SUBTRACT ARRCNT*250 FROM THE IRAWDAT INDEX IF THE INDEX 
BECOMES CREATER THAN ARRCNT*250. 

NO RM- LR EDPNT-ARRCNT*250*O VER FLW+S PNPRD 

IF (NORM .GT. ARRCNT*250) NORM-NORM- ARRCNT*250 

IDATA(LPNTCNT)-IBATA(LPNTCNT-l)-IRAWDAT(LREDPKI-ARRCtrr* 

+ 250*OVERFLW)+IRAWDAT(HORM) 


IF THERE ARE 250 ELEMENTS IN DATA, THEN WRITE OUT. 

IF (LPNTCNT .LT, 2 51) GOTO 700 
OUTPUT(l)-TIHE 

LOAD UP THE FIRST, THIRD, FOURTH, AND FIFTH CHANNELS OF THE 
DATA FRAME WITH JUNK. 

ALSO PUT IN AN OFFSET AS THE DATA OTHERWISE WILL BE CENTERED 
GROUND ZERO. 

OFFSET-0000 
DO 500 1-1,250 

IOUTPUT ( (1-1 )*5+l ) -IRAWDAT (I ) 

IOUTPOT ( ( 1-1 ) *5+2 ) -IDATA( 1+ 1 ) +OFFSET 
IOUTPUT (( I- 1 ) *5*3 ) -IRAWDAT ( 1*250 ) 
IOUTPUT((I-1)*5+4)»IRAWDAT(I+500) 

IOUTPUT ( (1-1 ) *5+5 ) -IRAWDAT( 1+ 7 50 ) 

500 CONTINUE 

IDATA ( 1 ) -IDATA( 2 51 ) 

WRITE(4) OUTPUT 
LPNTCNT-1 
TIME-TIME+. 050000 
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c 

C WHEN RKDPNT IS PAST A SUBBLOCK, THAT SUBBLOCK HAY 

BE DISCARDED ANDREPLACED BY A NEW ONE. THIS IS DONE 
VIA THE CETMEXT SUBROUTINE, 

’ 700 CONTINUE 

DO 800 I-l.AHRCHT 

IF ( (LREDFHT-ARRCNT*250*OVERPLW) ,EQ. 250*1) CALL 
CETNEXT(I,IRAWDAT) 

800 CONTINUE 

LREDPNT-LREDPNT+1 

LPNTCNT-LFNTCNT*1 

RETURN 

END 
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